THE UNIVERSIT Y OF MI CH IGAN COLLEGE OF ENGINEERING Department of Electrical Engineering SYSTEMS ENGINEERING LABORATORY Technical Report No. 2 RQA-1-THE RECURSIVE QUEUE ANALYZER V. L. Wallace R. S. Rosenberg ORA Project 07742 under contract with: ROME AIR DEVELOPMENT CENTER RESEARCH AND TECHNOLOGY DIVISION AIR FORCE SYSTEMS COMMAND CONTRACT NO. AF 30(602)-3953 GRIFFISS AIR FORCE BASE, NEW YORK administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR February 1966

TABLE OF CONTENTS Page LIST OF TABLES LIST OF FIGURES ABSTRACT 1. INTRODUCTION 2. NUMERICAL ANALYSIS FOR RQA-1 v vii ix 1 2.1 2.2 2.3 2.4 2.5 2.6 2.7 Conditions for Convergence Convergence for Continuous-Time Markov Chains Convergence for Discrete-Time Markov Chains Unique Convergence Estimation of Error Representation of Matrices Timing 4 4 5 6 6 7 8 9 3. EXAMPLES OF TRANSITION TABLES 3.1 Simple Server System 3.2 Fixed-Size Memory-Shared Computer Model ("FIXMEM' Model) 3.3 Random-Size Memory-Shared Computer Model ("RANMEM' Model) 3.4 Summary 12 13 18 23 28 4. PROGRAM DESCRIPTION 4.1 The Discipline Subroutine (DISCPL) 4.2 Stochastic Conversion Subroutine (STOCAL) 4.3 The Test Subroutines(TEST and TOUT) 4.4 The Estimation Subroutine(ESTIM) 4.5 The Iteration Subroutine, (TER) 4.6 The Iteration Strategy (MAIN and CVGTST) 4.7 Miscellaneous System Routines 4.8 Conclusion 29 29 29 31 31 33 34 55 36 5. PROCEDURES FOR THE USE OF THE PROGRAM 5.1 5.2 5.3 5.4 5.5 5.6 The Discipline Subroutine, (DISCPL) Program Control Parameters Structural Parameter Discipline Parameters Transition Table Parameter Initial Iterate 37 37 44 45 46 47 47 iii

TABLE OF CONTENTS (Concluded) Page 5J7 Format of Input Data4 5.8 Normal Output 48 5o9 Diagnostics of Abnormal Output 52 5.10 Format of Program Deck 53 -5.t-i.:ormat of the Punched Final Iterate 53 6. CONCLUSION 55 6.1 Timing Experience 55 6.2 Final Remarks 58 6 5 Acknowledgments 58 REFERENCES 59 APPENDIX A: LISTING OF RQA-1 60 APPENDIX B: SOME SUBROUTIT\ES FOR USE ING GENEERATING QUADRUPLES IN RQA-1 80 B1l Purpose 80 Bo2 Calling Sequences (MAD) 80 B 3 Restrictions 8i B 4 Discipline Program Structure 81 B.5 Example 82 iv

LIST OF TABLES Table Page Io Transition Table-"Simple Server" Model 17 II Internal-State Designation -"'EI.XMEM'1 Model 18 III. "FIXMEM' Model Quadruples 21 IV. Internal State Designation-"RANMENM' Model 23 V, "RANMEM' Model Quadruples 26 VI. "Simple Server" Discipline Subroutine 9 VII. "FIXMEM' Discipline Subroutine 40 VIII. "RANMEM' Discipline Subroutine 41 IX. Program Control Data 44 X. The Functions of PCH, PRSW, INIT, and TSW 45 XI. Structural Parameter 46 XII. Typical Output of RQA-l 49 XIII. "FIXMEM' Model56 XIV "RANMEM' Model 57 v

LIST OF FIGURES Figure Page 1. Fixed-size memory-shared computer model FIXIMEM). 19 2. Random-size memory-shared computer model (RANMEM). 24 3. Flow diagram for RQA-1. 30 4. DISCPL subroutine. 43 vii

ABSTRACT The Recursive Queue Analyzer, RQA-1, is a computer program designed to evaluate the equilibrium joint probability distributions of queue lengths and system conditions in very large, finite Markovian queueing systems. It will rapidly treat discrete-time and continuous-time systems having as many as 5000 states through numerical solution of equations of balance. The algorithm is generally very much faster and more accurate than equivalent Monte-Carlo simulation. This report is both a description of the organization of the program and a manual for program use. ix

1. INTRODUCTION The Recursive Queue Analyzer, RQA-1, is a computer program designed to evaluate the equilibrium joint probability distributions of queue lengths and system conditions in very large, finite Markovian queueing systems. This report is both a description of the organization of the program and a manual for program use. In the study of networks of queues and servers, the Markovian models represent an exceedingly important class of models. This program has been designed to facilitate the analysis of both discrete-time and continuous-time Markov chains having as many as 5,000 states. The primary design goal has been to provide a computation fast enough to encourage experimentation with the models in the study of system design. This has been achieved through efficient use of available (32K) high-speed storage in the computer and through careful program design. The program was written in the MAD language1 with selective use of the UMAP assembly language. Although this work has been motivated by the need for it in the design of large computer systems, the potential area of application is quite broad. For a continuous-time Markov chain with a finite state space, one can always write the Kolmogoroff differential equation At = At (1) Pt = Pt Q where p is a vector whose i-th element is the probability that the system is in the i-th state at time t, p represents the time derivative of p, and Q is a matrix of constants called the transition intensity matrix of the chain. The matrix Q is assumed to have the properties that ij = 0 (2a) for all states i, and qij > 0 (2b) for all i f j. An equilibrium probability distribution is, by definition, a solution to the equation 1

p Q = 0 However, any solution to Eq.. (3) is also a solution to p = p(AQ + I) (3) where A is a scalar constant and I is the identity matrix, RQA-1 is to solve Eq. (4). The purpose of For a discrete-time Markov chain with a finite state space, one can always write the Chapman-Kolmogoroff equation Pn+l = nA (5) where Pn is the vector whose i-th element is the probability that the system is in the i-th state at time n, and A is a matrix of constants called the transition matrix of the chain. The matrix A is assumed to have the properties that J a.. = 1, 1 zj (6a) for all states i, and 1- >ai (6b) for all i and jo A real matrix satisfying these properties is called a stochastic matrix. An equilibrium probability distribution in this case is a solution to the equation p = pA (7) Because of the obvious similarity between Eqs. (4) and (7), RQA-1 can conveniently so.3ve either equation through a suitable specification of the input. The most common application, however, is expected to be the continuous-time case, and, hence, a major part of the theoretical discussion and examples is devoted to that case. 2

The state space will usually originate as a multidimensional vector space with components corresponding to lengths of queues, numbers of occupied servers of a particular type, and similar quantities. However, since the number of such states is assumed finite, it is always possible to represent the state as a single variable i whose values are in correspondence to distinct values of the vector. The variable i is the state variable described in the preceding paragraphso Nevertheless, the elements pi of the solution to Eqs. (4) and (7) can be interpreted as joint probabilities of the variables represented in the original multidimensional state space, and appropriate marginal and conditional distributions can be calculated from them as desiredo This report will first treat, in Section 2, the numerical methods which are applied to the solution of Eqs. (3) and (7). This will be followed in Section 3, by a treatment of the technique used for the storage of the transition (intensity) matrix information in the form of what we call a "transition table" and several illustrations of the technique needed for the modeling of systems and the determination of transition tables Section 4 and Appendix A constitute a description of the structure of the program and the function of its component routines, while Section 5 explains in detail how the program is applied. A person whose only purpose in using RQA-1 is to obtain a set of equilibrium joint probability distributions should be able to skip Sections 2 and 4 without major loss, while a person desiring to alter or to transliterate the program for other computers or purposes will, of course, be most interested in those sections. The conclusion, Section 6, will discuss the performance of the program and some of the problems associated with its use.

2. NUMERICAL ANALYSIS FOR RQA-1 The RQA-1 employs an iterative procedure to determine the equilibrium probabilities for homogeneous continuous-time and discrete-time Markov chains, An equilibrium probability is found even when the chain is cyclic or when the chain does not have a unique equilibrium probability distribution. This section is concerned with the conditions for convergence and for uniqueness of solutions. It details the mathematical justification for the procedures used. 2o1 CONDITIONS FOR CONVERGENCE If G is an arbitrary m-th degree complex matrix, we will say that the iteration process p(n+l) p(n) G, (8) where p(n) is an m-dimensional row vector, converges if a nonzero vector JTm (n) Lim p(n) _ (9) n-oo existso We will say that the iteration process converges uniquely if it converges and if the vector, r, is unique (iae., independent of the choice of p(O) the initial iterate) Clearly, convergence of the process to - implies that T= rGo A set of sufficient conditions for the iteration process (Eqo (8)) to converge is that there be an eigenvalue of G at unity, that all other eigenvalues of G which are not equal to unity have modulus less than unity, and that all the unity eigenvalues of G have elementary divisors2 of first degree0 For, from Eqo (8) we may clearly write that p(n) = p(0) Gn (10) Then from the confluent form of Sylvester's Theorem,3 we may also write (making use of the above conditions) that 4

p(O) Gn =.. + P(O) Gn Xl(al Y1 + a2 Y22+ +o ar Yr ) + 0(\l (11) where Xi represents the i-th eigenvalue (assumed ordered so that = 2 = 02 0 = - r > Ihr+ll > ooo > liml ), r is the multiplicity of the dominant eigenvalue k1^ the ai are scalar constants, and the yi are m-dimensional rowvectorso If X1 = 1, Eq0 (11) evidently converges to (al Yl + ooo + ar Yr), a row-vectoro Note also that the rate of convergence is determined by kr+l, the subdominant eigenvalue 202 CONVERGENCE FOR CONTINUOUS-TIME MARKOV CHAINS For the continuous-time Markov chains, we are concerned with the convergence of the process p(n+l) = p(n) (AQ + I) (12) where Q is a transition intensity matrix. Let R = max lqiijo Then, if A < 1/R, the matrix (AQ + I) is found to be i a stochastic matrix, and, hence, it can be shown4 that every elementary divisor corresponding to a unity eigenvalue is of the first degreeo Since the degree of these elementary divisors is not dependent upon the value of A, this implies 1 that they are of first degree regardless of A as long as A < R Moreover, it is a direct consequence of a theorem of Gershgorin5 that the eigenvalues of a transition intensity matrix Q will lie on the disc of complex values y such that!7 - RI < R (13) where R = m x Iqii o The eigenvalues of AQ + I must, therefore, lie on a disc of complex values such that x- (AR + 1) | < A R. (14) Thus, if A is chosen so that A < 1/R, then there can be no eigenvalues of (AQ + I) having unit modulus which are not of unity value, and the iteration process (Eqo (12)) is convergent to a solution of Eqo (4)o 5

On the other hand, while a choice of A < 1/R is sufficient for convergence, it is not always necessary. Indeed, a choice of A slightly greater than 1/R can often lead to faster convergence than one slightly less, but one cannot know by any simple test when this is so. Hence, in RQA-1, A is chosen equal to.99/R to guarantee convergence for continuous-time modelso 2,3 CONVERGENCE FOR DISCRETE-TIME MARKOV CHAINS For discrete-time Markov chains, we are concerned with the convergence of p(n+l) p(n) A (15 where A is a stochastic matrix. As a result of the theorems cited above, A is known to have only elementary divisors of first degree associated with its unity eigenvalues and to have eigenvalues lying on the unit disc. In order to guarantE convergence, one may transform Eqo (15) into the iteration process p(n+l) = p(n) (99A +.01I), (16) noting that p = p(o99A + oOlI) (17) has the same set of solutions as Eq. (7), As a result, the eigenvalues of ( 99A + 01I) lie on the disc of complex values X' such that 1'. -.011 <.99 (18) which does not admit values of unit modulus not equal to unity. The iteration process (Eq. (16)) is used in RQA-1 for discrete-time models. 2.4 UNIQUE CONVERGENCE Under certain circumstances, a Markov chain may not have a unique equilibrium probability distribution. These are the cases in which the limiting probability distribution (lim pt or lim Pn) is dependent upon the initial distribution, p0o It is in ttioe circu1stances that the RQA-1 algorithm will not con 6

verge uniquely It converges, rather, to an equilibrium distribution which is consistent with the initial iterate, choseno This should cause no difficulty since such cases are rare in practice, and the user of RQA-1 will usually be aware of the dependence upon initial conditions through knowledge of his modelO In general, a sufficient condition for the convergent iteration processes above to converge uniquely is that it be possible (with nonzero probability) to reach each state from every other state in a finite timeo A transition matrix having this property (known as irreducibility) is called a primitive matrix Necessary conditions for uniqueness are somewhat strongero If nonuniqueness of a solution is suspected, a simple technique can be used to test the modelo The solution can be found by starting with an initial iterate having unit value for some (arbitrary) state and zero for all others. If the resulting solution has no states with a calculated zero probability, the solution is unique. If it has states with a calculated zero probability, then the solution process should be repeated using as an initial iterate a probability distribution which is nonzero wherever the previous solution was zero (within computational error) and zero wherever the previous solution was nonzero0 If the new solution is within the calculational error of the previous solution, it is the unique solution. If not, there is no unique solution. For more detail on the conditions leading to nonunique equilibrium distributions, Parzen, Gantmacher,4 or Chung7 are useful references0 2,5 ESTIMATION OF ERROR It can be shown that, if the elementary divisors of the subdominant eigenvalues Xr+l of a stochastic matrix G are of a degree m which is small compared to n, and if n is sufficiently large, then the n-th iterate can be approximated roughly by p(n) Gn p(O). + r+ n x (1 T( +X% n x (19) where x is a constant vector and it is the convergent solution0 Also, it can be shown that the difference between successive iterates can be approximated roughly by (n+l) (n) n - x (20) P fn+l) - pt uXs r+l(xr+l - l)n (20) A A fairly crude, but useful, estimate Xr+l of ['Xr+ll can be obtained from 7

1 1 m-1 1 A n n n n ^ A p(n+l) p(n)r+ll r+ 11 hr+1 ^ In - r-nl h~r+l - i I1n x n I.(21) since, for n sufficiently large, all but the first multiplicative term (kr+l) approach unity. Using this estimate, it can be argued from Eq. (20) that (n+l) (n) P -P n m-1 p -p+ — - A nml x (22) Ar 1 - Or+l which is an approximation to the error term in Eq. (19). Provided t at the +r+l calculated is less than unity and greater than about.5, (kr+1 - 1)1/ can be assumed to be near enough to unity, and the estimate so provided is,at the very least, a more realistic estimate of error than (p(n+l) - p(n)). 2.6 REPRESENTATION OF MATRICES In order to store information completely describing matrices A or Q, with A or Q having the requisite high degree, a special approach is required. Moreover, because it is desirable to have the matrix elements expressed as algebraic functions of a set of parameters which may assume different values at different executions of the program, this approach must allow expression-of the matrix elements in a literal form in the program. More will be said about the latter requirement in Section 35. For now, we will consider the problem of compressing the information needed to describe A or Q in such a way that computations are not appreciably slowed by the decoding required. Clearly, a 5,000-degree matrix, when stored as a two-index array, requires 25,000,000 locations of storage which is unreasonable for a "fast" program. However, both A and Q are generally sparse matrices (have mostly zero-valued elements) and will usually have a high degree of repetition of equal element values. Hence, a scheme of storage which lists location information along with value information is a plausible starting point. We construct a set of four vectors, together called a transition table, which implicitly define the matrix. Let us call them a,, 7, and B* and denote their i-th elements by ai, ~i, 7i, and Bi, respectively. The quadruple (1i ~i' 7i Bi) specifies one or more elements of a matrix in the following manner. *In the program listing, these correspond to the variables named "TRANS," "FIRST," "DEST," and "LAST," respectively. 8

The value of the element is ai, and its matrix co-ordinates are (pi, Yi)O Due to the repetitive nature of these matrices, the element ai may occur in other locations of the matrix with co-ordinates (pi + rS, y7 + rb) where 5 is a constant (fixed throughout the transition table) and r takes values 0, 1, 2, oo,(Bi - pi)/o In other words, the quadruple (ai, pi, 7i, Bi) specifies the occurrence in the matrix of an element with value ai at co-ordinates (Pi, 7i), (Pi + 5, 7i + 6),,oo,(Bi, 7i + (Bi - Pi))o This set of coordinates will be called the range of the quadruple. The set (ip, hi + 6, ~~, Bi) will be called the p-range of the quadruple The constant 5 will be called the period of the table0 The choice of a value for 5 will be discussed in detail in Section 3, Also in Section 3, we will carry out the construction of transition tables for three systems of varying degrees of complexityo These systems will be characterized as consisting of a single dominant queueing facility and a complex service facility, The degree of.the Q-matrix for such a system will depend on the maximum number of jobs allowed in the queue, This number will be a parameter in the transition table, and, in this way, a transition table of fixed size can represent a Q matrix of arbitrary ordero 2o7 TIMING Each iteration (Eqso (12) or (16)) is carried out in the following sequence~ n+1 (1) The vector to contain p is initially zeroedo (n) (n l). (2) al is multiplied by p( and accumulated into p (3) C1a is multiplied by p() and accumulated into P(n1+ if Pi + 6 < Bi (4) Step (3) is repeated until the range of the quadruple is exhausted (5) Steps (2) through (4) are repeated using ai, Ai 7i, B for i = 2,'3 ooo until all quadruples have been treated As a result of this procedure, the number of multiplications per iteration is equal to the number of nonzero elements in the transition matrix, In general, if the average number of different events which can occur at any time is E and the number of states is S, then the-number"of multiplications per iteration is ES. The convergence-error after n'iterates is O(Xr+l) Thus, if E is the desired error, the number of iterations needed will be I " logE/log|lr+l|, and the total number of multiplications required in the main loop is 9

M = ES loge/logl|r+lI (23) In contrast, a rule of thumb estimate for the standard deviation of a measurement of a probability p in a continuous-time Markov chain using simulation is s 2 (24) 7jT where T is the length of the simulation (in system-time) and y is the nonzero eigenvalue of the matrix Q having smallest modulus. For e to be equal to two standard deviations, this mean's that T- should be T = -, (25) and the number of random numbers generated would need to be K = 8 vp (26) L17E2 where v is the mean rate of occurrence of events~ It can be shown that IZI/v is usually of the order of -log|lr+ll which, in turn, is normally much smaller than unity. Hence, a relatively crude measure of the computation time of simulation and RQA-1 calculation is obtained from Eqso (23) and (26) by approximating the ratio K/M by K = 1 (27) KM 2( 27 ) M \e2 loge ES/ 2 Since e loge is very small for E <.01, e log e.1 43,5.01 2,170.001 145,000 10

we conclude that RQA-1 has a large advantage over simulation if desired accuracy is moderately high, This conclusion is strengthened considerably by the fact that a much greater amount of "Toverhead" calculation is required in simulations than with the RQA-1 algorithm. 11

30 EXAMPLES OF TRANSITION TABLES The preparation of the transition table is the most technical aspect of RQA-1 use. As stated in the previous section, the aim is to code a representation for a transition matrix which has the flexibility to include a wide class of problems and which makes efficient use of storage. The general format of the transition table has already been discussedo The purpose of this section is to illustrate the procedure by which the set of quadruples can be developed when a statement of the control rules and the structure of a Markovian queueing system is given. The examples are obviously not exhaustive, but it is hoped that they are sufficiently representative and that a user of RQA-1 familiar with them will have a minimum of trouble deriving transition tables for a wide variety of problems. Future reports will expand upon the process of modeling general, large queueing systems by Markov chains and by RQA-1 transition tables The procedure illustrated here is capable of being mechanized by computer programs to some degree once the principles are understood. In general, however, each queueing model represents a distinct programming problem, and it appears difficult to generalize to automatic techniques. For this reason, the present report will limit itself to essentials. More sophisticated general methods are the object of current research by the Laboratory. All three examples discussed in this section represent continuous-time Markov chains containing one main queue of jobs to be served by a service facilityo In each succeeding example, the service facility has an increasingly complex structure. The derived transition table in each case represents a transition intensity matrix Q& The pattern to be followed in transcribing a Markovian model for RQA-1 is: (1) Determine the state space and a suitable mapping of the states (which are almost always vectors) to a set of integers. (2) Determine the most useful value of 50 This will generally be the value of b which allows the transitions to be described with the least number of quadrupleso The choice of state-space mapping strongly influences what is possible here, and steps (1) and (2) may have to be attempted several times for a complex model. Computation time does not depend upon the number of quadruples used. Thus, the primary necessity here is to remain within available storage where space for only 2,000 quadruples is allotted. A secondary necessity is to keep the transition table preparation from becoming unnecessarily unwiedly. (3) Choose a sequence for generating the quadruples. This sequence must insure that every possible event.f- '.the model is considered for every 12

possible state of the model and that no event is considered more than once for any state, It must also insure that, when an event-state pair is considered, it is incorporated by extending the range of an already prepared quadruple whenever possible rather than creating a new quadruple, 3 1 SIMPLE SERVER SYSTEM The first and simplest example is illustrated in the diagram below: X ----- /^ \ --- j-Server N Jobs arrive as a Poisson process with mean rate k, The duration of job processing time is independent of arrival and is exponentially distributed with mean time 1/I,o Jobs are executed on a first-come first-served basis, If more than N jobs are present in the queue, arrivals are simply rejected, The first step is to define the state space for the system, Because of the simplicity of the example, the definition of state is artificial in this instance but provides a pattern for future use. The system state is defined by the number of jobs in the waiting line and by whether or not the service facility is occupied, Thus, a state is designated by the ordered pair (i,j) where i takes values in the set I = (0,1,2,, ooN) and j takes values in the set J = (0,1). In this case, the cardinality C of the set J is 2. It is apparent that not all ordered pairs (i,j) in the Cartesian product I x J represent states, In fact, the set [(i 0)| i >. 03 is not a set of states, since, physically, a queue can exist only when there is a job in service, In general, the Cartesian product I x J can be partitioned into two setst the set of states of the system, and the set of ordered pairs which are not states which we will call points, Thus, for this example, the set of states of the system is ((0,0), (0,1), (1,1), (2,1), o., (Nl)), and the set of points is ((1i,0), (2,O),, o o og(NgO) As stated previously, the state description within the program is onedimensional. A natural mapping, T, from the two-dimensional description (itj) to the mingle dimension, k, is k = T(ij) = Ci + j + 1 o (28) 13

Thus, the one-dimensiona- state description k runs from k = 1 corresponding to (i,j) = (0,0) to k = C(N + 1) corresponding to (i,j) = (N,C - 1) which is equal to (N.1), in this caseo In this example, the set of states is (12,4,6,, oo6,2N + 2}o At any time, the system can change state due to a job arrival or a service completion or can remain in the same state if neither of these occurso The transition table can be partitioned into three classes of quadruples: those associated with job arrivals, those associated with service completions, and those associated with no change, i.e., (1) Job arrivals, (2) service completions, and (3) no change. We can subdivide each of the above classes by the first dimension (a) i = 0; (b) i= 1; (c) i = 2, etC., and further subdivide them by the second dimension i) j =0 (ii) j = o A subclass is interpreted as a collection of transitions from a particular state caused by a particular event. For example, subclass (1) (b) (ii) represents the transitions from state (1,1) caused by job arrivals. One can then create quadruples of the form (a,7,yB) by sequencing through each subclass, making certain that the transition rate associated with each subclass corresponds to the a of a quadruple, and that the state associated with the subclass is in the p-range of the same quadruple (see Section 204). The state corresponding to the destination of the transition must be the y corresponding to that po The following detailed construction should clarify the above discussion. In this example, 6, the period of the transition table, is chosen to be equal to 2, since, as will be seen, the states to which equal rates apply differ most often by 2. It will usually be the case, when the model contains a main queue, that 5 is most conveniently chosen to be equal to C, the cardinality of the set J of all possible combinations Of values of.-al state variables other than the queue lengtho However, in general, 6 may be any other value (constant throughout the transition table) which allows the matrix to be described with a smaller number of quadrupleso 14

Beginning with the subclass (1) (a) (i), the first two parts of a quadruple can immediately be written, i 1 o The fact that a = X is obvious since class (1) is concerned with arrivals and jobs arrive at the system with a rate X.o = 1 is the result of applying the mapping T to the state- designated by (a) (i), namely (0,0), To determine y, we note that an arrival in state (0,0) produces a transition to state (0,1) and, therefore, 7 = T(0,l) = 2~ To determine B, assume B = r5 + P for r O 0; then if X applies for the transition 1 +2 it should next apply for the transition ioe, 35 4. But in the listing of states, 3 does not occur and hence B = P = l, Thus the first quadruple (x,1,2,1) applies only for the transition from state 1 to state 2 and, therefore, refers only to state 1, or subclass (1) (a) (i)o The next subclass is (1) (a) (ii)o We can immediately write a = X P = T(0,l) = 2 o An arrival in state (0,1) produces a transition to state (l14) and, therefore, y = T(ll) = 4o Again, assume B = rb + P for r f 0; then \ will apply for the transition 2 + 5 4 + i e,, 4 - 6 To determine B, it is necessary to find the maximum r such that r b + 7 < 5(N + 1) and such that the corresponding transition is a correct description of an arrival ioeo, in this case, r = N + 1 - 4/2 x N -,l In general, the maximum r is 15

desired for which the transition intensity a for states P + rb to 7 + rb is still valid. For all the examples to be treated, this will correspond to the maximum r for which y + rb is in the state space Therefore, B = 2(N-1)+2 = 2No Thus, the second quadruple (,2,4,2N) has the range. ((2,4),(4,6),ooo,(2N22N+2)) and, therefore, refers to the subclasses (1) (a) (ii), (1) (b) (ii), (1) (c) (ii' etco Since (b) (i), (c) (i), etco, are not states, all the subclasses of class (1) have been referred too A similar procedure now begins for class (2)0 The P for the quadruple associated with subdivision (2) (a) (i) would be 1, but, since there can be no service completions in the 1 state, the corresponding a is zero and requires no quadruple We next turn to (2) (a) (ii) for which = 2 and a = x, the service rate0 A service completion in state 2 produces a transition to state 1, ioe., 7 = lo Assuming B = rb + P, for r f 0, the next transition should be 4 o But this is impossible since a completion does not take state 4 into state 3 (which in fact, is not a state) Hence,the resulting quadruple is (p.,2,1,2) and refers only to subclass (2) (a) (ii). Since (b) (i), (c) (i), etco, are not states, the next subclass is (2) (b) (ii) for which a = pi and 3 = 4o The corresponding y is 20 B will be determined from the maximum r such that rb + 7 < 5 (N + 1) i.eo, r = N + 1 - 2/2 = N Therefore, B = r6 + B = 2N + 20 The quadruple ( p,4,2,2N + 2) has the range (4,2),(6,4), 00 o o(2N + 2, 2N)) which refers to subclasses (2) (b) (ii), (2) (c) (ii), etc,, and thereby completes class (2) The quadruples for class (3) can easily be constructed by using the property of Q that the row sums are zero, For all these quadruples, it is also the case that P = y, For class (3) (a) (i), it can be readily determined that a = -\ X = 1 7=1 since the only way to leave state 1 is through an arrival with transition intensity?o Thus, the practical rule for determing the a's of quadruples for 16

class (3) is that a equals minus the sum'of all transition intensities associated with leaving that particular stateo For the above quadruple, it is apparent that B = 1 since state 1 is the' only state for which there will be no,change if there are no arrivalso The next subclass is (3) (a) (ii) for which we write a = -(K + p) = - 2 7 = 2 o The fact that a = -(X + p) is an illustration of the above rule where the sum of the transition intensities is X +.o This quadruple will refer to all states which have the same properties as state 2, ioeo, states 4,6,ooo,2N, Hence, B is 2N, and the quadruple is (-( + i), 2,2,2N) For the final state, 2N + 2, the only transition intensity associated with leaving is p; therefore, the necessary quadruple is (-,2.2N + 2, 2N + 2, 2N + 2) This completes class (3), thereby completing references to all subclasseso The transition table is listed in its entirety in Table Io We will now go on to consider a more complex problem in considerably less detail. TABLE I TRANSITION TABLE — SIMPLE SERVER" MODEL 7 B Job- X 1 2 1 Arrivals X 2 4 2N Service 2 1 2 Completions 4 2 2N + 2 No - 1 1 1 Changes -x-2 2 22N -p,2N + 2 2N + 2 2N + 2 17

3 2 FIXED-SIZE MEMORY-SHARED COMPUTER MODEL ("FIXMEM' MODEL) This model was discussed in the paper by Fife and Rosenberg8 It is de picted in Fig. lo Jobs arrive at the loading queue as a Poisson process with mean rate ho Each job finds some number n of jobs in queue ahead of it. The queue is cleared on a first-come first-served basis, and the job enters the first stage of "service" where the mean rate of loading of jobs is il and the number of jobs being loaded is nl(n1 = 0 or 1)o If more than N jobs are in queue, the arriving job is rejected, After loading is completed, a job enters the second state of service, the processing phase. There are 5 t 'channels" in this stage, and n2 is the number which are occupied The mean rate at which jobs are processed in a single block or channel is pi2 A job leaves the system when its processing time ends. The principal complexity of the model is introduced by the following constraint: nl + n2 < 5 This implies that loading for no job can be begin while all five blocks are processing jobso The state description consists of the ordered triple (n,nl,n2) where n takes values 0,l,,o,,N, n1 takes values zero or 1, and n2 takes values 0,1,2,3,4,5, with the restriction that n1 + n2 < 5. Table II lists all the 11 possible occupancy configurations (nl,n2) and associates a value j to eacho TABLE II INTERNAL-STATE DESIGNATION-" FIXMEM' MODEL Internal nl n2 State j 0 0 1 0 1 2 0 2 3 0 3 4 o 4 5 0 5 6 1 0 7 1 1 8 1 2 9 1 3 10 1 4 11 18

JOB ARRIVALS AT RATE X / JOB COMPLETIONS n QUEUE nI LOADING PROCESSING Fig. 1. Fixed-size memory-shared computer model. 19

The values of j are here called "internal stateso" Thus, for this example, we have the set I = (1,2,ooo,N] and the set J = (1,2, oo,11) and, again, we take =6 C = 11 where C is the cardinality of Jo The mapping T for this example reduces to k = T(i,j) = 11 i + j The set of "points" is [(ij)li > O, j 1,2,354,5 ) Table. III lists the 47 quadruples which are necessary to refer to all the subclasses of this system. There are more subclasses for the example than for the simple server system since there are 11 internal states rather than 2o We will not carry out the complete derivation of this table but, rather, we will consider a few illustrative cases. We will derive the quadruples which refer to state (0,6) for all three classes For class (1), we can write a a 6 An arrival at state 6 will remain in the queue resulting in state (l,6), ioeo, state 17 0 Thus y = 17o It is obvious that later arrivals will also remain in the queueo Thus, we are looking for a B = w6 + B such that wS + y < &(N+l) where w is the largest integer for which this inequality holds. The above yields w < N + 1 - 17/11 Therefore, w = N - 1 This results in B = ll(N-1) + 6o The quadruple (A,6,17,11(N-1) + 6), which is number (6) in the listing, has the range ((6,17),(17,28),,,(11(N-1) + 6,11N + 6))o For class (2), we write a = 512 X - 6 03 6 since state 6 represents the configuration in which all 5 processing blocks are occupied, and jobs are processed at a mean rate of t1 per block, the rate for 20

TABLE III "FIXME4' MODEL QUADRUPLES No. a y B (a) Job Arrivals 1 1 7 1 2 2 8 2 3 3 ~ 9 3 4 4 10 4 5 5 11 5 6 6 17 ll(N-1) + 6 7 7 18 11(N-) + 7 8 8 19 11(N-l) + 8 9 9 20 11(N-1) + 9 10 10 21 11(N-1) + 10 11 11 22 11(N-1) + 11 (b) Service Completions 12 p2 2 1 2 13 212 3 2 3 14 32 4 3 15 42 5 4 5 16 52 6 5 6 17 52 17 11 11N + 6 18 k1 7 2 7 19 1 18 8 11N + 7 20N + 8 21 P1 8 3 8 22 Pl 19 9 iiN + 8 25 2p2 9 8 11N + 9 24 9 4 9 25 P1 20 10 11N + 9 26 352 10 9 11N + 10 27 P1 10 5 10 28 21 11 11N + 10 29 4 2 11 10 11N + 11 30 1 11 6 11N + 11 (c) No Changes 31 -x1 1 1 32 -X-p2 2 2 2 33 -X-2P2 3 3 3 34 -k -3 2 4 4 4 55 -x-4P2 5 5 5 36 -X-5p2 6 6 11(N-1) + 6 37 -5L2 11N + 6 11N + 6 11N + 6 38 -k-(1 7 7 11(N-1) + 7 39 -P1 11N +7 11N + 711N + 7 40 -\-pl-P2 8 8 11(N-1) + 8 41i lN + 8 llN + 8 11N + 8 42 -X-p1-2p2 9 9 11(N-1) + 9 43 -~1-2L2 11N + 9 11N + 9 11N + 9 44 — P1-35P2 10 10 ll(N-l) + 10 45 -ii-35 11pN + 10 llN + 10 l1N + 10 46 — l-, -4L2 11 11 11(N-1) + 11 47-Pp 11N + 1- 2 11 11 1N + 11 11N + 11 21

I blocks in use, where I~ 5 is 4t1. in this instance, I = 5, and, therefore, a = 5p2. A single service completion will result in a state in which I = 4, and, therefore, y = 5. If there is a queue, then a service completion in the processing block will result in the job at the head of the queue moving into the loading phase. Thus, two quadruples, namely (52.,6,5,6) and (5k2,17,11,11l+6), will be necessary to form the range ((6,5),(17,11),(28,22),...,(11N + 6, 11N)). Class (3) will require two quadruples for the same range: for the first a = -( + 5p2), and for the second a = -5p2. These quadruples are numbered (36) and (37), respectively. For the states 6,17,28,...,11(N-1) + 6, transitions are possible in two ways: an arrival with transition intensity X, and a service completion with transition intensity 5.2' For state llN + 6, a transition is only possible through a service completion. For a final case, we will derive the quadruples which refer to state (0,9) for all three classes. State 9 is the service configuration consisting of 2 jobs being processed and 1 job being loaded. Thus, there are two ways for service completions to occur. The quadruple associated with class (1) is number (9) in the listing, and, at this point, its derivation should be straightforward. For class (2), the relevant quadruples are numbered (23),(24) and (25). The a for quadruple (23) is 2p2 which is the transition intensity associated with a service completion in the processing blocks. The range of this quadruple is ((9,8),(20,19),...,(11N + 9, llN + 8)}. The service configuration for state 8 is shown in Table II. There are two quadruples necessary to describe service completions in the loading block. When there are no jobs in the queue, a service completion in the loading block results in an additional processing block being employed, i.e., state 4 results. When there are jobs in the queue, the job at the head of the queue moves into the recently vacated loading block and state 10 results. The class (3) quadruples are numbered (42) and (43). The range of quadruple (42) is ((9,9), (20,20),...,(11(N-1) + 9, ll(N-l) + 9)), and its a is -(k + 1 + 2p2), as expected. Quadruple (42) refers to state 11N + 9 with an a of -(1i + 2k2). These few cases should be sufficient to indicate the structure of the entire transition table. The next example will introduce an additional complexity in that service completions will depend more closely upon the nature of job arrivals. 22

3.3 RANDOM-SIZE MEMORY-SHARED COMPUTER MODEL ("RANMEM" MODEL)* The model depicted in Fig. 2 bears a resemblance to the previous model discussed, but it is only a superficial one. As before, the queue will be cleared on a first-come first-served basis. There will be three types of jobs arriving, distinguishable by the number of service blocks or storage units they require. The jobs will arrive as a Poisson process with mean arrival rates X\, \2, X3 for one, two, and three storage unit jobs, respectively. These rates are independent of time and state of the system. For convenience, we set X = 1l + +2 + 3. We assume that the system is nonpreemptive, i.e., a job once begun is completed. Each type of job will have a basic mean service rate when it is being processed alone. These service rates are pl, k2', p3 for one, two, and three storage unit jobs, respectively. However, when two jobs are being processed simultaneously, the mean service rate at which each job is completed is assumed to be one-half the basic mean service rate. For example, if a oneunit and a two-unit job are being processed, the mean service rate for the one-unit job is (1/2)k1 and for the two-unit job (1/2)X2. Similarly, for three one-unit jobs, each will have a mean service rate of (1/3)k11. Upon completion of any type of job, all blocks being used for that job become free simultaneously. This model is particularly interesting because it requires the evaluation of conditional probabilities (conditional upon the state of the system) in order to determine the correct transition intensities. This is true because jobs having different properties are allowed to be mixed in the queue, and the nature of the first job in the queue is dependent statistically upon the condition of the system. The state description consists of the ordered quadruple (n,nl,n2,n3) where n is the number of jobs in the waiting line, and nl,n2,n3 are the numbers of one-, two-, and three-unit jobs in service, respectively. Table IV lists all TABLE IV INTERNAL STATE DESIGNATION-" RANMEM' MODEL nl n2 n3 Internal State, j 0 0 0 1 0 0 1 2 0 1 0 3 1 1 o 4 1 O 0 5 2 0 0 6 3 0 0 7 *This model corresponds to one suggested by Dr. R. V. Evans, now of UCLA, for multiprogrammed time-shared computer operation. 23

| __ I JOB ARRIVALS AT RATES (XI X2, X3) (- I 1 F1 t~.= S WL IVLWl - I ] I ] JOB COMPLETJONS SERVICE AT RATES COMPLETIONS (,41,,2,7 3 ) Fig. 2. Random-size memory-shared computer model (RANMEM) 24

the possible occupancy configurations of which there are 7. For this example, the set J is (1,2,...,7) with 5 = C = 7, and the set I is (0,1,...,N). The mapping T reduces to k = T(i,j) = 7 i + j. We will again call the value of j the "internal state" of the system. The set of "points" is ((i,j)li > O,j = 1}. Table V lists the 60 quadruples which are necessary to refer to all the subclasses of the state space. We will consider only one illustrative case. We will derive the quadruples which refer to state (n,3),n = 0,1,...,N, for all three classes. State 3 represents a two-unit job in service. The arrival of a one-unit job with transition intensity Xl produces a transition to state 4. The arrivals of either a two-unit or a three-unit job with intensity k2 + X3 do not change the internal state but do increase the queue length. Both these quadruples (5) and (6) refer only to the single state 3, since for state 10, the arrival of any type of job (with intensity X = \1 + k2 + k3) produces a transition to state 17. Then quadruple (14) will have the range ((10,17),(17,24),...,(7(N-1) + 3, 7N + 3)). Thus, three quadruples are necessary for the class (1), subclasses (n,3), n = 0,1,...,(N). In state 3, a service completion will produce a transition to state 1 with intensity 2.' Quadruple (18) will have only the single element (3,1) in the range, because the transitions resulting from service completion of a twounit job are quite different for state 10. In state 10, the single job in queue can only be a two-unit or a three-unit job (a one-unit job cannot remain in the queue when there is space for it in service. If, for example, only a two-unit job in service, the one-unit job cannot remain in the queue). We have here the situation in which the identity of the- first job in the queue depends partially on the type of job in service. For example, if there is a two-unit job in service, the first job in the queue cannot be a one-unit job. If there is a one-unit job in service, the first in queue cannot be a one-unit or a two-unit job but must be a three-unit job. This dependency will introduce complexities into the form of the transition intensities. If there is a two-unit job in the queue, a service completion will result in a transition with intensity t2X2/(X2 + X3); similarly, for a three-unit job in the queue, the transition intensity will be >2X3/(k2 + k3). In the first case, the \2/(x2 + %3) factor is the probability of a two-unit job being in queue given that there is a two-unit or a three-unit job in queue. Quadruples (27) and (28) account for these cases. Quadruple (2) has a complete range because, given a three-unit job at the head of the queue, the remaining composition of the queue is unimportant. But quadruple (27) has only the single element (10,3) in its range, since transitions from state 17 depend on the nature of the first two jobs in the queue if the first in queue is a two-unit job. The only two cases to consider for state 17 are that the first in queue is a two-unit job and the second in queue is not a one-unit job, and that the 25

TABLE V "RANMEN" MODEL QUADRUPLES No. a 7 B 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 k1 k2 X3 X1 xi k2 + k3 X2 x3 ~2 + X3 k2 k (a) Job Arrivals 1 1 1 2 3 3 4 5 5 5 6 6 7 10 12 13 5 3 2 9 4 10 11 6 4 12 7 13 14 17 19 20 1 1 1 7(N-1) + 2 3 3 7(N-1) + 4 5 5 5 6 6 7(N-1) + 7 7(N-1) + 3 7(N-1) + 5 7(N-1) + 6 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 13 112 1l/2 P2/2 I11 3x3 /x 1132/ ^3x1/% I 2x2/( 2 + x3) 52x3/(2 + + ) tl(X3 + 2)/ 2)3 /2) (41ll + 112~2)/2% 112X1/2% 1l5x3/(X2. + 35) '1x2/(x2 + x3) 1i ILl(x2 + 3)/x P3x2( 2 +2%3) /x 2p3x23T/x 2L3X-/k2x2/X 12x2xi/(X( 2+x3 ) 1,2xx2 + \%)/32 -x -2 -11/2-.2/2 -L1 -1l -k -1 -111 -k-|-l (b) Service Completions 2 3 4 4 5 '6 7 9 9 9 10 10 11 11 11 11 13 13 12 14 14 16 16 16 16 17 17 18 18 23 23 (c) No Changes 1 7N +2 2 7N + 3 3 7N +4 4 7N + 5 5 7N +6 6 7N + 7 7 1 1 3 5 1 5 6 2 3 5 3 2 10 12 4 6 12 4 2 13 7 10 4 12 16 10 4 13 7 13 7 2 3 4 4 5 6 7 7N + 2 9 9 10 7N + 3 7N + 4 7N + 4 7N +4 11 7N + 6 7N +6 7N + 5 7N + 7 7N + 7 7N + 2 7N +2 7N + 2 16 7N + 3 7N +3 7N + 4 7N + 4 7N + 2 7N + 2 48 49 50o 51 52 53 54 55 56 57 58 59 60 1 7N + 2 2 7N +3 3 7N +4 4 7N+ 5 5 7N +6 6 7N+ 7 7 1 7N + 2 7(N-1) + 2 7N +3 7(N-l) + 3 7N +4 7(N-l) + 4 7N +5 7(N-l) + 5 7N +6 7(N-1) + 6 7N + 7 7(N-l) + 7 26

second in queue is a one-unit job. In the first case, the transition intensity will be 42x2/1 = 2 2 +3 X i.e., k2 multiplied by the probability that the first job in queue is a twounit job and the second is either a two-unit or a three-unit job. The transition will be made to state 10. The complete quadruple is number (42). In the second case, the transition intensity will be 2 1 -2 2 + X3 ' i.e., p2 multiplied by the probability that the first job in queue is a twounit job and the second is a one-unit job. Thus, after the service completion two jobs move into service resulting in state 4. Number (43) in Table Vb is the complete quadruple. For this problem, it will be advisable to determine the class (3) quadruples by examiniation of those classes (1) and (2). The sum of the transition intensities referring to state 3 is M2 + %l + k2 + x3 resulting in the transition intensity -(K + p2). To determine the range of the quadruple with this a, we will consider a's for states 10 and 17. For state 10, the following is the sum of the transition intensities. x + 22/(x2+ x3) + 2%3(x2 + x) This also reduces to X + kt2. For state 17, the following is the sum of the transition intensities: X + k2%3/(X2+ ) + 2X22/X + 2X) 2X1/(X(X2 + 3) ' This again reduces to X + p2. Thus, quadruple (52) will refer to all states except the final state. For the final state, the a will be simply -p2 since quadruple (14) with a = k does not refer to the final state. A few additional observations may be useful. In quadruples (19) and (20), the transition intensities p1 and.2 are each divided by 2. This is a result of an assumption made earlier in which it was stated that, if two jobs were run simultaneously, the mean rate at which each would be completed would be 1/2 the basic mean rate. When two like jobs (i.e., one-unit jobs) are being run, the mean rate for each is (1/2)[1i but the mean rate for a service completion will be the sum, namely 1.. Quadruple (39) has a factor 2 associated with the transition intensity. This arises from the fact that two like transition intensities have been added in order to arrive at this result. 27

3 4 SUMMARY It should be noted that all three transition tables have a common featureo Each represents a system with the size of the state space limited by the maximum number of allowed jobs in the queue. As we stated previously, this is only one kind of system for which a transition table may be usefully employed. Additional complexities can be introduced by having elements of the quadruples expressed as functions of some or all of the parameters of the queueing modelo As illustrated, the method for deriving a transition table provides an orderly procedure for determining a complete set of quadruples. What is actually done is to enumerate all events which can occur for each state. The division of events between arrivals and service completions in the examples is artificial and done only for convenience in ordering tasks. For problems in which some events cannot be identified as either an arrival or a service completion, some other division will be suitableo In any case, one will continue to enumerate the transitions for each event and each state. Of course, sequences other than those suggested in the examples could be followed in systematic evaluation of the quadruples. In some instances, the construction of a partial state diagram will greatly aid the process of determining the transition table. For complex systems, the state diagram may prove very difficult to draw and may not provide much insight into the regularities of the state space. Once again, we stress the point that examples have been presented to suggest possible forms for transition tables. For most problems, it will be the user s responsibility to devise transition tables with those properties most useful to his special purposes.

4. PROGRAM DESCRIPTION The purpose of this section is to discuss programming details of the RQA-1 program. For the most part, the block diagram (Fig. 3) and the MAD listings (Appendix A) of the program provide full documentation of the program. However, a few remarks will be made here to aid in their interpretation. We shall begin these remarks by discussing the roles of the major subroutines. More detail of the procedures to be followed by a user, of definitions of user specified parameters, and of diagnostic printouts is provided in Section 5 which acts as a sort of "user's manual" for RQA-1. 4.1 THE DISCIPLINE SUBROUTINE (DISCPL) The function of the Discipline Subroutine is to prepare the transition table which describes the transition probabilities or intensities of the problem to be solved. This subroutine is unique to each application or model solved by RQA-1, and hence, must be written by the user. The form of the transition table internal to RQA-1 is a set of four vectors called TRANS, FIRST, DEST, and LAST corresponding to the a, P, y, and B respectively, of the quadruples discussed in Section 3. The i-th element of each of the four vectors together form the i-th quadruple. The DISCPL subroutine is basically a program for the evaluation ofthese vectors in terms of parameters which may be altered at execution time. These parameters are the counterparts of the X, pt, N, etc., of the examples in Section 3 and appear only in DISCPL. In addition to preparing the transition table, DISCPL has the responsibility of reading all parameters supplied by the user at execution and of providing essential control information descriptive of the queueing model being analyzed. Sample DISCPL routines for the three models used as examples in Section 3 can be found in Section 5 (Table VI, VII, and VIII). 4.2 STOCHASTIC CONVERSION SUBROUTINE (STOCAL) The role of this subroutine is to determine whether the transition table in memory represents a transition intensity matrix Q (hence, a continuoustime model) or a stochastic matrix A (hence, a discrete-time model). If it 29

MAIN DISCPL STOCAL ENTRY DETERMINE APPROPRIATE TRANSFORMATION OF MATRIX TRANSFORM MATRIX TO CONVERGENT STOCHASTIC MATRIX CRETURN TEST. ENTRY TEST FOR ERROR IN FIRST OR LAST TEST IF STATE ROW SUMS-I & POINT ROW SUMS=0 TEST IF POINT COLUMN SUMS=0 TESTS YES SUCCESSFUL RETURN NO DIAGNOSTICS (TOUT.) TOUT. ENTRY PRINT DIAGNOSTIC AIDS RETURN ESTIM. ENTRY CALCULATE INITIAL ITERATE PERFORM I ITERATION (ITER.) RETURN ITER. ENTRY P-AP RETURN PROUT ENTRY PRINT PROBABILITES AND MARGINAL PROBABILITIES RETURN Fig. 3. Flow diagram for RQA-1. 30

is the former, STOCAL converts the transition table to represent the stochastic matrix( (99/R)Q + I), where R = max qii. If it is the latter, STOCAL converts the transition table to represent the new stochastic matrix (.99A +.01I). These conversions assure the convergence of the resulting iteration processes to an equilibrium distribution of the specified model (cf., Sections 2.2 and 2.3). The determination of the type of matrix represented is accomplished by determination of the element of the TRANS vector having least value. If it is equal to or greater than zero, the transition table represents a stochastic matrix. If, on the other hand, it is negative, the transition table represents a transition intensity matrix, and the value of the element having least value is equal to -R. 4.3 THE TEST SUBROUTINES TEST AND TOUT The purpose of TEST is to determine whether the current transition table does indeed represent a stochastic matrix, and whether it adheres to the necessary properties of a valid transition table. Using the transition table, the elements corresponding to each row of the matrix are summed and successively stored in a prezeroed vector. If the matrix is stochastic, the rows corresponding to states will have a one in the associated position of the vector and these corresponding to points will have a zero, In addition, all elements in columns corresponding to points must be zero. In the TEST s.ubroutine, it is determined whether or not the above conditions exists. If they do not exist, the subroutine TOUT is called, and various diagnostics are printed to help determine the error(s) in the transition table. For each quadruple in a valid transition table, the following equation must holds B - = rb for r e (0,1,2,...) (29) If this condition is not met for all the quadruples of the transition table, the subroutine TOUT is called, and appropriate diagnostics are printed. 4 4 THE ESTIMATION SUBROUTINE ESTIM Before beginning the iterations, an initial iterate must be selected. The RQA-1 main program allows this iterate to be selected in four different ways by the user: it can be read in as input data either in binary or floating point form, it can be taken to be the final iterate of the previous run (which has remained in storage), or it can be calculated by the subroutine called ESTIM.

As its name suggests, ESTIM is a subroutine for estimating a satisfactory initial iterate and is based on a generalization of the simple server system discussed in detail in Section 3.1. For a simple server model, the probability of queue length i is proportional to pi, where p = x/p. (e.g., ref. 10)o ESTIM prepares an initial iterate so that the k-th element of the iterate is proportional to (p*)k, where p* is a constant. If G = (gij) is the matrix represented by the transition table, then p* is calculated as 1 P* = (-) (30) where = gi (31) i,j i> j i,j i>j li - Ji ij: giji 0~ v = (33) and T is the total number of nonzero, off-diagonal elements of the matrix. For the simple queue, these. expressions reduce to p* = kX/. The only real justification for applying this calculation to more complex models is that it has been found to be quite effective, particularly when'.the model being solved originally had a multivariate state space and the mapping to the linear state space had a major queue of the system as the most significant variable in a Cartesian mapping. ESTIM calls the iteration subroutine once before returning to the calling program, producing a vector having nonzero values only in states. The vector supplied to the calling program will not have a unit sum. Hence, subsequent normalization by the calling program is required. 32

4.5 THE ITERATION SUBROUTINE, ITER Whereas the main program and all the other subroutines are written in the procedure-oriented language, MAD, the ITER subroutine is written in the machine-oriented language, UMAP.* The purpose of this is to make this subroutine as efficient as possible. The only argument of this subroutine is the number of iterations to be carried out. Thus, a typical call of the subroutine would be EXECUTE ITER.(N) where N is the number of iterations. One iteration consits of the multiplication of the current iterate by the matrix. The subroutine requires two vectors in erasable store: one for the last iterate, and one for the iterate being formed. A switch is set which indicates which is the current iterate, thus removing the necessity for continually transposing newly computed iterates. Because of the importance of this subroutine, we would like to discuss its operation in more detail. The two vectors mentioned above are actually stored in one vector which is represented in the program by the symbol V. This vector has 10,000 locations, thereby allotting 5,000 to each vector. To indicate which half of V contains the current iterate and which half will be used to store the newly computed iterate, two symbols, VFLAG and CVFLAG, are defined in the program. These symbols can only assume the values 0 or 1 and are constrained to have complementary values, i.e., when VFLAG = 1, CVFLAG = 0. The current iterate is always referred to by calling V(VFLAG,I). When VFLAG = 0, the first half of V contains the current iterate. When VFLAG = 1, the second half of V contains the current iterate. To illustrate how a typical iteration would proceed, assume that VFLAG = 1. This means that the current iterate occupies the second half of V. During the process of multiplying this vector by the matrix, values of the newly calculated iterate are stored in the first half of V which had been zeroed previously. When this single iteration is complete, VFLAG is set equal to 0, and the second half of V is zeroed. Thus, the current iterate occupies the first half of V and will be correctly referenced by calling V(VFLAG,I). As mentioned above, the multiplication process is carried out by sequencing through the transition table. For a particular quadruple (ai,yi,Yi,Bi), ai is multiplied by the value of the current iterate at location Bi' and the result is added to the value at location 7i of the iterate being formed. This process is continued for the range of the quadruple, i.e., until ai is multiplied by the value of the current iterate at location Bi with the result being added to the value at location 7i + Bi - Pi of the iterate being For users who do not have a UMAP assembler available, a MAD version of ITER is also included in the listings of Appendix A. Its use in place of the UMAP version will approximately quadruple the execution time of RQA-1.

formed. The program then proceeds to (ai+lPi+l,7i+lBi+l) and repeats the process until all quadruples are exhausted. This type of vector-matrix multiplication is quite efficient because only nonzero elements of the matrix are referenced. Also, no significant penalty in computing time results from the manner of describing the matrix. 4.6 THE ITERATION STRATEGY (MAIN AND CVGTST) Upon completion of the calculation and normalization of the initial iterate, the main program produces a first call to the convergence testing subroutine CVGTST. This first call of CVGTST simply sets the underflow traps so that underflows produce normal zeros, and returns to the calling program (MAIN) with an index corresponding to a nonzero element of the initial iterate. The main program continues to iterate until successive values of the specified element have an absolute difference which is less than the convergence criterion E, or until the specified maximum number of iterations has been taken. At that point CVGTST is called again. If we define p(n) to be the n-th iterate and pi(n) to be the value of the i-th element o the n-th iterate, the maximum difference between successive iterates e n) is given by e(n) = max I p(n) (n-l) max |p - Pi The iteration process is said to have converged after n iterations for error criterion E f e(n) < e. For every call of CVGTST after the first, CVGTST calculates et) and tests for convergence. If the process has not converged, this subroutine returns to the main program with the integer iM corresponding to the state for which the absolute value between successive iterates is largest. The idea is not to calculate iM for each iteration but to use an initially calculated iM as a coarse test and to recalculate only when this coarse test is satisfied. Then the main program iterates until the absolute difference of the iM-th element is less than e or the specified number of iterations has been taken. Another call of CVGTST then results, and the process is repeated. On the other hand, whenever the tests within CVGTST are satisfied, CVGTST prints appropriate comments and error estimates, resets the underflow error trap, calls the output subroutine OUT2, and transfers out (via an "error return") to the beginning of the main program for a calculation of another model (if data are available). 34

The purpose of this two-level testing procedure is to avoid the expense of calculating iM for each iteration by using a value of iM in a coarse test for convergence and recalculating iM (using CVGTST) only when the coarse test has been satisfied. 4.7 MISCELLANEOUS SYSTEM ROUTINES For the sake of completeness, the library routines explicitly called by RQA-1 will be listed here with a brief mention of their function. (a) BNBCD.(X) converts its argument X from binary to binary coded decimal (Hollerith) form suitable for punching or printing characters. (b) BPUNCH.( ) punches the values of the list fin column binary format.* (c) BREAD.(/ ) reads the values of the list from column binary cards. (d) ELOG.(X) computes the natural logarithm of X. (e) ERROR. produces a call to the monitor system with a dump of memory if requested. (f) FTRAP. sets the underflow trapping procedure so that underflows are replaced by normal zeros until NTRAP is executed. (g) NSDBCD.(X,Y) provides for punching the BCD values X,Y in columns 72-80 of the next card punched. (Sequencing of the last characters of this field is accomplished by The University of Michigan Executive System.) (h) NTRAP. resets the underflow trapping procedure to its normal condition so that an underflow produces an exit to the monitor system (as in ERROR). (i) SYSTEM. produces a call to the monitor system. (j) TODAY.(DATE,YEAR) gives to DATE.and YEAR BCD values corresponding to the day, month, and year that the program is executed. (k) ZERO.(gL) sprays normal zeros into all locations specified by the list J. Actually, in The University of Michigan system, it writes binary card images on the "punch output" tape, for offline punching. 35

4o8 CONCLUSION We would like to make some remarks on the program as a whole Except for the two subroutines DISCPL and PROUT,* the program operates on a state space which is one-dimensional. Thus part of the program needs no information with respect to the structure of the states, i.e., their multidimensional nature. The transition table created by DISCPL describes a onedimensional process, but DISCPL also contains some information about the structure of the states which is passed on to PROUT. For example, if the state designation were two-dimensional, DISCPL would contain values for the upper bounds of each dimension. In the examples of Section 3, the upper bound for the internal state variable is C which is fixed, and the upper bound for the queue state variable is N which is put in as data. The upper bounds are needeld only in the subroutine PROUT to prepare a suitable output formats PROUT was written specifically for a two-dimensional state description and presents the final iterate in terms of such a state description. Thus, the two subroutines DISCPL and PROUT can be considered usercontrolled. DISCPL must be supplied by the user, while PROUT may be replaced by the user if he desires output formats other than those built into RQA-1 or if he wishes to perform some additional caculation before printing results. *PROUT is the program,called by OUT2, which prepares and prints the printed state probabilities and marginal probabilities. 36

5 PROCEDURES FOR THE USE OF THE PROGRAM This section describes procedures for the use of. RQA-o1 It'is concerned primarily with the data which must be supplied by the user and with the interpretation of the output produced by the programo It is assumed here that the program is used with The University of Michigan Executive Systemo In normal use, the object deck of RQA-1 is followed by the DISCPL subroutine, a "$DATA" card, and the input data Since the DISCPL subroutine must be prepared by the user, we will give a detailed discussion of its preparation. The data are in five classes, divided into two distinct groups of cards. In order of their appearance following the $DATA card, these are GROUP 1 (1) Program control parameters (2) Structural parameters (3). Discipline parameters (4) Transition table parameter GROUP 2 (1) Initial Iterate. Following the description of the programming of a DISCPL suvbroutine will be a discussion of the input data, the form of both no.rmal and abnormal output, the structure of the program card deck, and the format of binary cards occurring as punched output 5ol THE DISCIPLINE SUBROUTINE, DISCPL The chief input to RQA-1 provided by the user is the discipline subroutine. This user-generated program provides RQA-1 with a complete description of the model to be solved in the form of code which generates the transition table, In addition to correctly describing his model, the user must insure, while writing this subroutine, that the number of quadruples generated does not exceed 2,000 and that the linear index of every state does not exceed 5,000O* The discipline subroutines must carry out certain necessary functions for RQA-1: * Failure to do the former will usually result in a violation of low-core storage protection in the University of Michigan system and probably worse difficulties in other systems. Failure to do the latter will simply result in an error return to system with suitable diagnostics printed out. 37

(1) Specifying the complete transition table (2) Reading or specifying the Group 1 parameters (3);,Specifying a name )for the model. The transition table consists of four vectors called TRANS, FIRST, DEST, and LAST corresponding respectively to the ca, P, y, and B of Section 3. The values of elements of these vectors will usually be dependent upon the values of discipline parameters which are supplied at execution time and, hence, cannot be treated as simple numerical vectors to be read as data. Rather, they are computed by the execution of program segments contained in the discipline subroutineo This facilitates examination of the behavior of a model as parameters are varied, since the program may be simply rerun with different numerical values of the parameters supplied. Normally, the discipline subroutine must contain all the statements within the dashed blocks in Fig. 40 The first block contains the instruction to read the Group 1 data, and in this section, it will be assumed that this is done by the simplified input/output statements in MAD. The second block defines the model name ("name" should be replaced by any desired twelve characters, including spaces) and sets up necessary communication to the remainder of RQA-1 through PROGRAM COMMON. Group 1 parameters which are not read as data may be declared or calculated in DISCPLo The Group 1 parameters will be defined later. The vectors TRANS, FIRST, DEST, and LAST of the transition table can be set up in several ways, only one of which is exhibited in the listings of the DISCPL subroutines shown in Tables VI, VII, and VIII for the examples of Section 5. This method involves the use of two list-manipulating statements in MAD, namely: SET LIST TO and SAVE DATA.1 SET LIST TO TRANS sets up a vector named TRANS, and SAVE DATA Ol,, O2,.oo, n causes al, a2.o ~ on to be the n entries of this vector where the a's can be any arithmetic expression acceptable to MAD. While the illustrated technique is quite useful, other methods will also be found suitable to particular problems. For example, the quadruples may be generated in a series of THROUGH loops when some systematic evaluation technique is available or the number of quadruples to be generated depends on data Appendix B contains a discussion of a group of subroutines, contributed by J. H. Jackson, which have been found to be useful in the generation of quadruples under quite general circumstances. On the whole, the degree of innovation which is possible is limited chiefly by the user's desires and needs. In determining the quadruples of the transition table, the period DELTA (6) can be treated as a constant or a parameter. If it is a constant, as is the case in the examples of Section 3, it is set in DISCPL, but, if it is to be treated as a parameter, it can be defined in DISCPL or read in as data.

TABLE VI "SIMPLE SERVER" DISCIPLINE SUBROUTINE SCOMPILE MAD D ISCQ0OO 000787 09/09/65 2 09 11.6 PM MAO (09 AUG 1965 VERSION) PROGRAM LISTING....... FXTERNAL FUNCTION DISCOOlO 001 ENTRY TO DISCPL. DISC0020. 002 INPUT PARAMETERS ARE L, MU, EPSI, NQ _ PRINT COMMENT $8 INPUT PARAMETER.....SCO03 *003. 1I S $ DISC0040 *003 READ AND PRINT DATA _______.._...._...__ DIS050_ 004 MAXJ = 2 DISC0062 *005 DELTA = 2 DISC0063 *006 SET LIST TO TRANS DISCO080 *007 S'_AVE DATA LL MUMU.MUL -L,-L-MU-,MU __...__._ _ ___ _.._ISCQ.90__.008 _ SET LIST TO FIRRT DISCQ100 *009 SAVE DATA 1,2,2,4,1,2,2*NQ+2 DISCOlO *010 SET LIST TO DEST DISC0120 *011 -_ SAVE DATIA-L2. 4, 1, 2, 2 22 NQ+2 __ _.. DI..SCQ13Q..*...... 2... 12. SET LIST TO LAST _..__ DISCO.140 *013 SAVE DATA 1,2*NQ,2,2*NQ+2,1,2*NQ,2*NQ+2 DISC0150 *014 VECTOR VALUES NAME = $SIMPLE MODELS DISC0160 *015 INTEGER FIRST,DELTLAST, INITPCHPRSWDELTA,MAXJ,MXITE,N DISC0200 *016 PROGRAM COMMON TRANS(2001),FIRST(2000) DEST(2000) LAST(20_00) _ DISC220 *017 1 PCHPRSW,MXITER,EPSIDELTAINIT,RUN,TSW,NAME(l),MAXJ DISC0230 *017 _FUNCTION RETURN _.....__ _ __. _ISC.0270 _.._ 18. END OF FUNCTION DISC0280 *019

TABLE VII "FIXMEM' DISCIPLINE SUBROUTINE SCOMPILE MAD,EXECUTE,DUMPI./O DIUP,PRINT OBJECT,PUNCH OBJECT OISCOOOO:-.~' ' ' 0)0791 09/09/65 4 40 21.8, PM MAD -039 'AU6,- 1965 VERSION) PROGRAM LISTINS 000 000 000, 0 EXTERNAL. FUNCTION ENTRY TO DISCPL. INPUT PARAMETERS ARE L, U1l, NU2,,EPSI, NQ PRINT COMMENT $8 INPUT PARAMETER 1 S $ READ AND PRINT DATA MAXJ = 11 _ __ ____. __ _ DELTA = 11 ""S fT O -TR-" A'S --- —-, ---- --------— ' --- -- - SET LIST TO TRANS SAVE DATA L, LtLL,LL, LtLiL,L, 1 NU2,2*NU2,3*~U2,4*MU2,5*NU2,5*NU2,NU1,NU1,NU2,NUI, 2 NU1,2*NU2NUNUlNUl,3*NU2,NUlMU1,4*NU2,NU1, _____ 3 -L,-L-NU2,-L-2*NU2,-L-3*NU2, 4 -L-4*NU2,-L-5*NU2,-5*NU2,-L-NU 1,-NU1,-L-NU1-NU2,-NU1-NU2, 5 -L-NUl-2*NU2, -NU1-2*NU2,-L-NU13*NU2,-NUl-3*NU2 -L-NUl-4*NU2, 6 -NU1-4*N2 SET LIST TO FIRST_________.. SAVE DATA 1,2,3,4, 5, 6,t 7 8,9, 1,11,2,3,4,5,6,17,7,18,8,8,19,9, 1 9t,23,10, 1,21,11,11,l,2,3,4,5,6, 1llNQ+6,7,11*NQ+7,8,11*NQ+8,9 2,11*NQ+9, 1,'11*NQ+10, 11,11 NO+ll SET LIST rTO DEST SAVE DATA 7,8,9,10,11,17,18,19,20,21,22,1,2',3,4,5,11,2,8,7,3, 1 9,8,4,13,9,5,11, 1,6,l,2,3,4,5,6,1llNQ+6,7,11*NQ+7,8,11*NQ+8, 2 9,11*NQ+9, 10,11*NQ+10,11, 11*N 1 _ _ __ _ _ _ __.SET LIST TO LAST SAVE DATA 1,'2,3,4, 5, 1.l Q-1)+6, 11* (NQ-1)+7, 11(NQ-1)+8, 1 11 ( Q-1 )+9, 11( N Q-1)+1i, 11._( Q-)+11,2,3.4,5,6, llNQ+6,7, __ 2 11*NQ+7,11*NQ+8,8,11*NQ+8,11*~d+9,9,11*NQ+9,11*NQ+10,10, 3 11*NQ+13,11*NQ+11, 11*NQ+11, 12,3,4,5,11*(NQ-1)+6,11*NQ+6,. 4 11*(NQ-1)+7,11*NQ+7,l11*(Q-1)+8,11*NQ+8,11*(NQ-1)+9,ll*NQ+9, 5 11*. (Q-1 )+10, ll*Q+l, ll*(NQ-1)+11, l*1NQ+ll1 _ VECTOR VALUES NAME = $FIXMEM I3DELS INTEGER FIRST,DEST,LAST, INIT,tPtHPSWDELTA,MAXJMXITER,NQ PROGRAM COMMON TRANS(2001),FIHST(2000),DEST(2000),LAST(2000), 1 PCH,PRSr,MXITER,EPSI,DE-TA, INIT,RUN,TSW,NAME( 1),MAXJ FUNCTION RETURN END OF FJNCTION DISCOOlO c001 Dl SC0010 *001 DISCO020 *002 DISC0030 *003 DISC0040 *003 DISC0050 *004 DISC0062 *005 DISC0063 *006 DISC0080 o007 DISC0090 *008 DISC0100 *008 DISCo110 *008 DISC0120 *008 DISC0130 *008 OISCO140 *008 DISC0150 *008 DISC0160 *009 DISC0170 *010 DISC0180 *010 DISC0190 *010 DISCO200 *011 DISC0210 *012 DISC0220.012 DISC0230 *012 DI SC0240 *013 DISC0250 *014 DISC0260 *014 DISC0270 *014 OISC0280 *014 DISC0290 *014 DISC0300 *014 DISC0310 *015 DISC0350 *016 DISC0370 *017 DISC0380 *017 DISC0420 *018 DISC0430 *019

TABLE VIII "RANME"' DISCIPLINE SUBROUTINE $COMPILE MAD DISCOOO0 001195 09/13/65 12 58 6.7 PM ' MAD —WT-9 AUG-'AT-9'65 VERS'TON )' PRRO GRAM LISTING...... EXTERNAL FUNCTION ENTRY TO 'DISCPL. -N-INPUT PARAMETERS ARE L1, L2, L3, L, NU1, NU2, NU3, EPSI, NQ PRINT COMMENT $8 INPUT PARAMETER 1 S $ READ AND PRINT DATA MAXJ = 7 DELTA = 7 DI SCOOOO CI SCOO10 DISC0020 DI SC0040 DI SC0040 DISCC062 DI SC0063 *001 *C02 *003 -*003 *C04 *005 *006 -i=F-j SET LIST TO TRANS SAVE DATA L1,L2,L3,L,L 1,L2+L3,L,Ll,L2,L3,L,L2+L3,L,L,L,L, 1 NU3,NU2,NL1/2.,KU2/2.,NU,NUI, 2 NU1,NU3*L1/L,NU3*L2/L,NU3*L1/L,NU2*L2/(L2+L3),NU2*L3/(L 32+L3),NU1*(.L3+L2)/( 2.*L),NU2*L3/(2.*L),NL1*L1/(2.*L)+NU2*L2/( 4 2.*L),NU2*L1/(2.*L),,NUl*L3/(L2+L3),NU1*L2/(L2+L3),NT'I,NU1*(L2 5 +L3)/L,NU1*L1/LNU3*L2*( L2+L3)/(L*L), 2.*NU3*L2*L1/(L*L),NU3*L 6 1*L3/(L*L),NU3*L1*L1/(L*L),NU2*L2/L,NU2*L2*L1/((L2+L3)*L),NU2 7 *L1*(L2+L3)/(2.*L*L),NU2*L1*Ll/(2.*L*L),NU3*L1*L1*(L3+L2)/(L* 8 L*L),NU3*L1*L1*L1/(L*L*L) SAVE DATA -L, 1 -NU3,-L-NU3,-NU2,-L-NU2,-NU1/.-NU2/2.,-L-NU1/2.-NU2/2,-NU1, 2 -L-NU1,-NU 1,-L-NU1,-NU 1,-L-NU1 DISC0070 DI SC0080 DISCO090 DISCO 100 DISCO110 DI SC0120 DISC0130 DISCO140 DISC0150 DISC0160 DISC0170 DI SC0180 DI SC0190 *007 *008 *008 -008 — *008 *008 *008 * 008 *009 *009 *009 SET LIST TO FIRST SAVE DATA 11,1i,2,3,3,4,5,5,5,6,6,7,1.',12,13,2,3,4,4,5,6,7,9, 1 9,9,1,,1::,11,Ll,11,11,,1,13,12,14,14,16,16,16,16,17,17,18,18, 2 23, 23, 1, 7! 3+2, 2, 7*N0+ 3, 3, 7*1NQ+4,4, 7*NQ+5,5,7*NQ+6,6,7*NQ+7,7 SET LIST rO CEST SAVE- CATA 5,3,2,9,4,1,1,1,6,4,12,7,13,14,17,19,2Ci.,1,,3,5,1,5 i,6,2,3,5, ',2,1C,12,4,6,12,4,2, 13,7,10,4,12,6,10,4,13,7,13,7,1 2, 7:)+2, 2, 7*f,)+,3, 3, 7*NQ+4, 4, 7*NQ+5, 5,7*NQ+6,6,7*NQ+7,7 DISCO200 DISC0210 DISC0220 DISC0230 DISC0240 CISC0250 DISC0260 DISC0280 *0 f-' -. *011' *011 *011 *0'2 - - *013 *013 *013

TABLE VIII (Concluded) SET LIST TO- LAST DISC029-0o *T4 SAVE DATA 1,1,1,7*(NQ-1)+2,3,3,7*(NQ-1)+4,5,5,5,6,6,7*(NQ-1)+ DISC0300 *015 1 7,7*(NQ-1)+3,7*(NQ-1)+5, 7*(N.Q-1)+6,2,3,4,4,5,6,7,7*(NQ )+2,9 DISC0310 *015 2,9,10O,7*(NQ )+3,7*(NQ )+4,7*(NQ )+4,7*(NQ )+4,11,7*(NQ ) DISC0330 *015 3 +6,7*(NQ )+6,7*(NQ )+5,7*(NQ )+7,7*(NQ )+7,7*(NC )+2, DISCC340 *015 4 7*(NQ )+2,7*(NQ )+2, 16,7*(NQ )+3,7*'NC )+3, 0ISC0345 *015 5 7*(NQ )+4,7*(NQ )+4,7*(NQ )+2,7*(NQ )+2 DISC0350 *T15 SAVE DATA 1, 7*N+2,+2,7*(NQ-1 )+27*NQ+3,7*(NQ-1)+3,7*NC+4,7*(NQ- DISC0351 *016 1 1)+4, 7*iN+5,7*(NQ3-1)+5,7*NQ+6,7*(NQ-1)+6,7*NQ+7,7* (N-1)+7 DISC0352 *016 VECTOR VALUES NAME = $RANMEM MODELS DISC0360 *017 INTEGER FIRST,DEST,LAST,INIT,PCH,PRSW,DELTA,MAXJ,MXITER,NQ OISCO400 *018 PROGRAM COJMMCN TRANS( 2jO1),FIRST(2'Ci ),DEST(2CC ),LAST(200), DISC0420 *019 1 PCH,PRSW,MXITER,EPSI,OELTA, INIT,RUN,TSW,NAME'(I),MAXJ DISC0430 *019 FUNCTION RETURN DISCC470 *020 ENC CF FUNCTION DISC0480 *021 ro

EXTERNAL FUNCTION ENTRY TO DISCPL PRINT COMMENT $8 READ AND PRINT DATA INPUT PARAMETERS $ Portion of program defining Group 1 parameters which are not read as data Portion of program defining transition table in'terms of DELTA, TRANS, FIRST, DEST, LAST O 0 0 0 O 0 0 a o 0 a 0 0 * 1 O0 0D a Q a 0 a 00 VECTOR VALUES NAME = $ "Name"$ INTEGER FIRST, DEST, LAST, INIT, PCH, PRSW, DELTA, MAXJ, MXITER PROGRAM COMMON TRANS(2001), FIRST(2000), DEST(2000) LT(2 ), PCH, PRSW, MXITER, EPSI, DELTA,.iNcIT, RUN, TSW, NAW (l), MAXJ FUNCTION RETURN END OF FUNCTION Figo 4 DISCPL subroutineo

We will complete this discussion with a brief reference to the three DISCPL subroutines listed in Tables VI, VII, and VIIIo In all these cases, the SET LIST TO and SAVE DATA statements have been used. All these examples have been formulated in terms of a single queue and a complex server~ The parameter NQ is used in the transition tables to represent the maximum queue lengtho Of course, any integer-valued variables introduced by the user must be declared. For these examples, DELTA was chosen equal to C (the cardinality of the set J) and so the following pair of statements is used: DELTA = x MAXJ = x where x is 2,7,11 for the SIMPLE QUEUE, RANMEM MODEL and FIXMEM MODEL respectively, (MAXJ is a structural parameter (Group 1) which, in these cases, is also equal to the cardinality of Jo) Group 1 parameters are defined in the next three subsectionso 5 2 PROGRAM CONTROL PARAMETERS The program control parameters are.named PCH, PRSW, MXITER, EPSI, INIT, and TSWo The mode, preset value, and range of permitted values of each of these variables are given in Table IX If these parameters are not specified on any particular run, they have, for that run, the value last supplied on a previous run. If no previous run has specified their value, they have the "preset" value shown TABLE IX PROGRAM CONTROL DATA Preset Other Name Mode Value Values PCH INTEGER 1 0,1 PRSW INTEGER 1 0,1 MXITER INTEGER 100 > 0, < 10 EPSI FLOATING 10-4 > O, < lo POINT INIT INTEGER 0 0,1,2,3 TSW INTEGER 1 0 1 44

EPSI is the value of the convergence test criterion, Thus, when the maximum absolute difference between corresponding elements of successive iterates is less than or equal to EPSI, the convergence test is satisfied~ MXITER is the maximum number of iterations to be allowed in the calculationo If this number of iterations has not been satisfied, the calculation is terminated and output supplied~ Table X lists the significance of each of the values of the variables PCH, PRSW, INIT, and TSWo TABLE X THE FUNCTIONS OF PCH, PRSW, INIT, and TSW Name Value Function PCH 1 - Final iterate is punched in binary formato O - Final iterate is not punchedo PRSW 1 - Final iterate is printed by execution of PROUTo O - Final iterate is not printedo PROUT is not calledo INIT 0 - The initial iterate is calculated by the subroutine ESTIMo 1 - The library subroutine BREAD is called and an initial iterate in binary format is read ino This input is the punched output of a previous run. 2 - The initial iterate is read in by READ DATA statement and is in floating point format 3 - The final iterate from the previous run is used as the initial iterate for the present runo TSW 1 - The transition table is tested to assure its stochastic character and other properties (TEST is called)o O - The transition table is not testedo 5o3 STRUCTURAL PARAMETER Structural parameters convey dimensionality information to be output printing subroutineso For the output printing subroutine provided (PROUT), there is only one such parameter, MAXJo This parameter indicates the maximum extent of the least significant coordinate J of the two-dimensional output array of proba-.45

bilities. In PROUT, a state with index K will be treated as a state with vector value (I,J) as determined by the mapping K = MAXJ * I + J + 1 corresponding to Eq. (28). Thus, MAXJ corresponds to the cardinality of the set (J] described in Section 3. The preset value and the range of MAXJ are given in Table XI. If the preset value is used (by not specifying MAXJ) or if a value greater than 2500 is specified, the printed output will supply the results ten-to-a-row and fiftyto-a-block. TABLE XI STRUCTURAL PARAMETER Preset Name Mode Value Restriction MAXJ INTEGER 10 MAXJ < 2500 5.4 DISCIPLINE PARAMETERS These data consist of the required parameters for the particular queue discipline being run. They are defined by DISCPL only. If they are not floating point constants, their mode must be declared. For the examples discussed in Section 3, such parameters are the L's and NU's. For the FIXMEM model (Table VII), a data card with discipline parameters would be L = 1.0, NU1 = 1.5, NU2 =.25,... In general, these parameters are the independent variables involved in the quadruples of the transition table. 46

5o5 TRANSITION TABLE PARAMETER This parameter is DELTA, which is an integer. Its use was discussed in Section 51o. It is treated essentially like a structural parameter, but, while structural parameters refer to the state-space structure, DELTA is ac.tually a part of the description of the transition table, It has a preset value of 1o 5,6 INITIAL ITERATE As seen in Table X, when INIT is 1 or 2, the initial iterate is read in as data. For INIT = 1, the initial iterate must be presented in column binary. format. It is usually the set of cards containing a final iterate which was punched in binary mode in a previous run as the result of setting PCH = 1, The form of this iterate is described in Section 5.11. The entire set of output cards as described there must be used as input for the initial iterate, i.e., the header card followed by the binary datao For INIT ' 2, the initial iterate must be supplied in floating point form to be read by a READ DATA type statemento Unspecified values are automatically zero. Before being used, the initial iterate will be normalized to a probability vector. Therefore, at least one element must be specified. Typical data cards would be V(l) =.25, V(3) =.50, V(10) = o25* or v(l) =.20, o15, o,.05, v(50) o20,.15, oO.,.o5* Values can be given for any V(I) where I runs from 1 to the largest integer denoting a state These data are Group 2 data and must appear after the Group 1 data of the run. 5.7 FORMAT OF INPUT DATA There is no specified order for the Group 1 data. These data are in the format of cards to be read by a READ DATA)statement in MAD and must be followed (according to that format) by "11"t 7

Following the Group 1 data will be the Group 2 data, if required (i.eo, if INIT = 1 or INIT = 2)o Following this will be the data for the next run, if further runs of the model with different data values are desired. Any number of sets of data may be stacked in this wayo The run number sequence which identifies all output of RQA-1 can be started (restarted) at a user-specified value on any run by insertion of a card just ahead of the Group 1 data for the run with the format THE INITIAL RUN NUMBER IS XXX o This statement must start in the first column of the card, and -words must be separated by a single spaceo The XXX may be any non-negative integer of three or less digitso If this is not done, runs will be numbered starting with 1 for the first run after loading. The following may be a typical set of data cards for one run: EPSI =.0005, MXITER = 200, INIT = 2 L = 1, MU1 = 1.5, MU2 = o25* V(l) =.5,.5,.1,.1* Alternatively, another set may be: THE INITIAL RUN NUMBER IS 75 INIT = 1 L = 1, MXITER = 50, MUl = 1.5, MU2 = o25* (Binary card deck as described in Section 5o11) a o. 9 o Note that, in this case, EPSI will be the same value as it was on the previous run (or it will be o0001 if this is the first set of data) and that, in both cases, PCH, PRSW, and TSW are the same as they were in previous runs (or their preset values if this is the first set of data) 5,8 NORMAL OUTPUT The normal output of RQA-1, illustrated in Table XII, should be selfexplanatory. However, it should be noted that the probabilities are printed in terms of a two-dimensional state designation where PROB(I,J) is the limiting state

TABLE XII TYPICAL OUTPUT OF RQA-1 RQA RLSULTS__ FOR FIXMEM MOQDEL- DI.S.C.IPLINE. SEPT 1, i965 VERSIOI JF RiQ-1.RUN.__...._EP.19.65. INPUT PARA4METERS L=1.,NUl=1.5,NU2=:. 5, NQ=15*...THE-_ [ILAL IITERATE.WAS T... ERM INED BY THE S.URDlUTINE._ESI M1.. ****************************** SUCCESSFUL RUN ****************************** CONVERGENCE TEST CRITERION (EPSI).CC:C?^ MAXIMUM ALLOWED ITERATIONS (MXITEER) = 1O ITERATIO'S TAKEL = 89 MAXIMUM DIFFERENCE 3ETWEEr4 LAST TWO ITERATES = 1. -E-34 ESTIMATED ORDER CIF ERO-IJR =.1E-.3 (APPR(XI'IlATE CONVERGENCE FACTOR =.9016573) THE SOLUTION VECTOR HAS BEEN PUNCHEli

TABLE XII (Continued) THE FOLLOWING ARE THE LIMIIING STATE PROBABILITIES_.PROLQBIJ I. PROB(I) IS THE SUM OF PROB(I,J) OVER ALL STATES WITH THE VALUE I FOR THE FIRST VARIABLE _ t PR II PRO I) I 1_PR(IR tJ- _ __ __ _ - - -_ = 0 1 2 3 4 5 6 7 8 9 0.52348530.04460908.08872072.08792652.05749297.02716861.00869784.02980799.05949457.05947242.03980388.02029071 _ _._ __ __ __ __ — - ____ _...... ----.. —.-. ---- -- --------- - - _. ---- 1.1526^A7Q.nQ000n000n.-nnnnnnn 00000000.00Q00000 -n0000 non -00.Rqi47.02000nn3s.- 04119641 -..040.695529 -..02792383.01496318 2.10608881.0000000.00000000.00000000.00000000.OCOOOOO0.00712090.01348743.02730851.02796113.01950914 -.01070170..___ —. __._ --—,._._ 3.073001n 0 _-00000000 -OOO0O000OO.0000QQ__- O000n000 0-0.00C00n0 0- 2 005951 0 —!70...01858519 011916128 -...1347094-.00744561 4.04970614.00000000.00000000.00000000.00000000.00000000.00364714.00617654.01261922.01303470.00916726 _, D D ~d6'.o 5 0 61 2. ___a.__........................................... _ _......._............. _......................... -........ -- -._5 __0345 35!t _. _ QO.OQO Q...Q...__G.iOC-.Q...C0OCC_ 024_74174.04..17181.G.U851817. 00877507.006142_41.00337126 6.02223266.0000000.O9000000.0000000.000G0000.OCQ00000.00164044.00280297.00569988.00583454.00405082 ____dl _ __.^00??0401 _0022.Q.._-1_ _.... 0 __.01457771 _.00000 _00.Q0 QQO..00000000. OOOCQO.OQQ.QQ 0..0106571. 00186844 3. 00377258.00382645.00262z885.00141568 8.00942488.00000000.00000000.00^0000.00000000.CCGCCOOO.00067974.00123296.00246581.00247318.00167895.......................... _ 8 4............ -9 00.Q600726 -....30300 o00.0000000.COCCG OO 0C000. CCCCO.0 'CGoQ.00042634.000804G6.0C1 158983 cDC157486. 00105586.00055633 10.00377696.0)000000.00000000. 0.300000n.03CGCC000.oCOCCCO.000263b0.00oC51756.0J310166.00098847.00065495......-.......00034182 11.00234679.00000300.00000000.00000000.COCCOO3.000CC.OO0.0Q 0116089.300032868.00063387.C0061297.i0040237.33020801 12.00144560.0000OOOO0.00000300.00o.COO.OCo CC00000 3 000 300009656.00020615.r.%39372.C*0037850.C 024682.30012387,13.00088689.300330o33.oo00o 00o.00000000 0.G000.OCCCCCO.oCiCGCooC.00005556.300312857.0024633.00023900.CC014833.03336910 14.00054715..O3000000 003J00030 OC.3-'0.0..OOOOO").OCCCO '.00029t4.00308311.00 17104.000,14868.;0008031.0C003438. 15..00033959.03000000 0 C000000 0.GOOCOOo3.OCGCOOcco.C COOO"J.i C 0)1674.03~.08897,. 0011723.00007530.G00003186. 0C 000949

TABLE XII (Concluded) ____ _'ROB(J) IS. THE SUM OF PRDAB(I,J) OVER ALL STATES W.. ITH THE VALUE J FOR THE SECOND VARIABLE CHECKSUM PROB(J) J =0 1 2 3 4 5 6 7 8 9.99999920.04460908.08872072.08792652.0"J5 749297.0271 6861.04C48632.9 ': 2 83293. 182823Or.-.182583.12699661. 0 6 7 6 9 6 6............................................................ \-n

probability of state (I,J). The marginal probabilities PROB(I) and PROB(J), which are the sums of PROB(I,J) over J and I respectively, are also listed. Thus, for a fixed value of I, PROB(I) is printed and followed on the same line by the values of PROB(I,J) as J runs over its allowed values. CHECKSUM, the sum of all the elements of the probability vector, may not be exactly equal to 1 due to roundoff. 5.9 DIAGNOSTICS OF ABNORMAL OUTPUT A run can terminate without supplying answers to the given problem for several reasons; there may be no input data, the number of allowable states or entries in the transition table may be exceeded, or there may be an error in the transition table. In each of these instances, appropriate statements an diagnostic aids are printed. This termination is caused by execution of the MAD instruction EXECUTE ERROR. which causes a return to the operating system proceded by a dump of storage if requested. The following is a further explanation of a few of these diagnostics which may not be self-explanatory: (1) THE FOLLOWING QUADRUPLES HAVE AN ERROR IN FIRST OR LAST. This is printed if LAST(I)-FIRST(I) is not an integer multiple of DELTA. This is inconsistent with the rules for forming a transition table. (2) THE FOLLOWING QUADRUPLES HAVE AN ERROR IN DEST. This is printed if there is a transition specified into a point, and there is no other transition with that point as a source. (3) THE FOLLOWING ROW SUMS, SUM(I), COF THE STOCHASTIC MATRIX ARE INCORRECT. The appearance of this statement indicates that the matrix represented by the transition table (at this point, it should be stochastic) is not stochastic. If a row sum corresponding to a point is nonzero, then a transition is being defined for something other than a state. If a row sum corresponding to a state is less than one, then a possible transition is being overlooked. If it is greater than one, an extraneous transition has been included. 52

5.10 FORMAT OF PROGRAM DECK Two formats will be described; one for which the DISCPL subroutine is to be compiled, and the other for which it has been compiled and, similar to the rest of the program, is now on binary cardso For the first case, the following is the deck format: 2 yellow ID cards "1$ COMPILE MAD, EXECUTE, (etc.)" card DISCPL MAD deck "$ BINARY" card RQA-1 binary deck "$DATA" card DATA For the second case, the following is the deck format: 2 yellow ID cards "$ BINARY, EXECUTE (etc,)" card RQA-1 binary deck including DISCPL binary deck "$DATA" card Data 5.11 FORMAT OF THE PUNCHED FINAL ITERATE As discussed in Section 5.2, when PCH = 1, the final iterate will be punched on binary cards. To illustrate the form of this output, we will assume that the name of the model is SIMPLE QUEUE, the date is MAY 25, 1965, and this is the 15th run. The first card punched, called the header card, will have the following form: SIMPLE QUEUE RUN 15 25 MAY 1965 SIMPLE Q. Since, in the University of Michigan system, binary cards have only columns 73 -80 interpreted (printed on the top of card), this card will be identified by the first 8 characters in the model name appearing in columns 73-80. The next card will contain the linear index corresponding to the maximum state, called SCOPE, punched in binary mode. Columns 73-80 of this card will.contain 25 MAY 01 ioeo, the day, the month, and the number 01o The succeeding cards will contain the values of the final iterate in binary mode. The number of such values will be SCOPEo Columns 73-75 of these cards will contain

015 (i.e., the run number) Columns 78-80 will contain the numbers 001, 002, etco (i.e., a sequential numbering of the cards). Each final iterate punched by RQA-1 will automatically have this output card format. 54

6. CONCLUSION The purpose of the RQA-1 program is to evaluate equilibrium joint-probability distributions of queue lengths and system conditions for very large, finite Markov queueing systems. Although only these probabilities are obtained by the program, it is possible to estimate from them a wide variety of other performance statistics. The expected waiting and throughput times are easily determined as well as the probability distributions of any of the state variables (e.go, queue lengths, channel occupancies) by taking appropriate sums of state probabilitieso Special programs can be written to determine the above statistics which will have as their input the probability distribution calculated by the RQA-1 program. As part of this section, we will present some timing data for computer runs of the queueing models discussed previouslyo 6.1 TIMING EXPERIENCE Results are shown in Table XIII(a) for the FIXMEM model with X and-. 1 fixed and p2 varied. Table XIII(b) gives the results for one run of the model with 611 states. It is interesting to note that for each model the number of iterations per second varies inversely as the number of stateso Table XIV presents some results for the RANMEM model for which the discipline parameters are fixed but EPSI and the number of states are varied. For the problem with 77 states and EPSI = 10-4, the number of iterations per second is 483,5 For the problem with 357 states and EPSI = 10-4 this number is 10 The ratio of states is 357/77 = 4.6, and the ratio of iterations per second is 48o3/10 = 4.8. It should be noted that these results have been achieved for very sparse matrices In fact, for the FIXMEM model, there is an average of four nonzero elements per row. Thus, for a 611 state model, the ratio of the number of nonzero elements to the total number of elements is (4 x 611)/(611 x 611) =,0065, Using the equation developed in Section 2, we can compare the number of multiplications required for a solution of a typical FIXMEM model by simulation Vith the number of multiplications required by the numerical methods employed in RQAA-l. or the FIXMEM model with t2 = 53333, the EPSI satisfied was 4058 x 10 which corresponds to an e, in the equation, of about 10-3 A simulation for a model with these parameters would require at least 100 times as many multiplications as the numerical procedure or about one minute as compared to one-half second. For the FIXMEM model with p2 =.5 and EPSI = o6429 x 10-, a simulation would require about 150 times as many multiplications or about 4-1/2 minutes as compared with 2 seconds for an E = 103lO 55

TABIE XIII "FIXMEf' MODEL (The initial vector is estimated for each run) =, X = o = 1.5 EPSI Number Iteration Number of Total Maximum 112 /_ 10-4) of of Iterations Time Probability ( 10-4) Iterations (Seconds) per Second (Seconds) (a) Number of States = 101 1.0 0.4884 73 2.4 30 4 356 0.1251 0.5 0.6429 64 1..8 3555 3.6 0.0801 0.3333 0.4058 23 0.6 38.4 2.4 000525 0.2857 0.9508 63 1.8 35.0 356 0.0411 TOTALS 223 6.6 33 o8 13 2 (b) Number of States = 611 0.5 1.0 54 8.4 6.4 17.4 0.0783

TABLE XIV "RANMEN' MODEL (The initial vector is estimated for each run) X1 = 1.5 2 = = ~3 X3 =.2; 1 = 2, g2 = 1.5, 13 = 1.2 Number 3E PSI: Number Iteration Number of Total Maximum of States x 10-4 of of Iterations Time of States x 10-4 Probability Iterations (Seconds) per Second (Seconds) 77 1.0 58 1.2 4803 3.6 0.3828 77 0.1 102 2.4 42.5 4.8 0.3844 77 0.0478 115 3.0 381.3 5 4 0.3845 TOTALS 275 6.6 41.7 13.8 357 1.0 60 6.0 10.0 1538 0.3812

6o2 FINAL REMARKS This work represents the first stage in exploiting numerical solution of stochastic models for the analysis of large systems. The major disadvantage of the technique is currently the relatively great difficulty of preparing correct discipline subroutines for complex models. It is hoped that this disadvantage can be overcome by the development of automatic programs to do this work from a very formalized description of the model to be solved (i.e., a translator for a problem-oriented language)o With that accomplished, we believe the RQA-1 and its successors represent a very powerful tool indeed for system analysis. The programs listed have used the May 1, 1965, version of MAD and the May 5, 1965, version of UMAP. These programs are known as the September 1, 1965, version of RQA-lo They are designed to operate with The University of Michigan Executive System, although an effort was made to make the program as independent as possible. We believe that with the information presented here a potential user could adapt the program to any comparable system, 6 3 A.KNOWLEDGMENT The authors are particularly indebted to Professor Ro V. Evans whose original suggestions prompted this work, to Mr- D. W. Fife who made many valuable suggestions on its application, to Dr. K. B. Irani for his critical examination of this manuscript, to J. H. Jackson for contributing the work in Appendix B, and to E. S. Walter and D. L. Mason for help to prepare portions of the. final RQA-1 program..58

REFERENCES 1. Arden, Bo, Bo Galler, and R, Graham, MAD-The Michigan Algorithmic Decoder, University of Michigan Press, Ann Arbor, 1964. 2. Gantmacher, F. Ro, Matrix Theory, Vol. 1, Chelsea, New York, 1959, po 142. 35 Frazer, R. Ao, W. J. Duncan, and A. R. Collar, Elementary Matrices, Cambridge, University Press, New York, 1960, p. 83. (It follows from the definition in 3 and 2 that a root of multiplicity r has degeneracy r if its elementary divisors are of first degree.) 4. Gantmacher, Ro F., Applications of the.Theory of Matrices, Interscience Publishers, New York, 1959, Theorem 10, p 102. 5. Varga, Ro So, Matrix Iterative Analysis, Prentice-Hall, New York, 1963, pp. 16-17. 6. Parzen, E., Stochastic Processes, Holden-Day, San Francisco, 1964, ppo 248-258, 7. Chung,., Ko L., Markov Chains with Stationary Transition Probabilities, Springer-Verlag, Hamburg, 1960. 8. Fife, Do, and R. Rosenberg, "Queueing in a Memory-Shared Computer," Proceedings of the 19th National Conference of the ACM, Philadelphia, 1964 9. Evans, Ro Vo, "Several Queueing Systems Suggested by Computer Organization Problems," Internal Memo. 10o Saaty, To Lo, Elements of Queueing Theory, McGraw-Hill, New York, 1961, pp, 38-400

APPENDIX A LISTING OF RQA-1 This program was written in the MAD and- UMAP languages and processed in The University of Michigan Executive System. The listings are complete for the September 1, 1965, version of the RQA-1 program. The listings for the DISCPL subroutines of the examples of Section 3 are given in Section 5o Two versions of ITER are supplied, one in UMAP and the other in MAD. 60

$GOMPILE MAD EXECUTE,.DUMP I/O DUMP.PRINT OBJECTPUNCH OBJECT MAINCOOO 002300 09/20/65 10 43 25.0 PM NAD (09 AUG 1965 VERSION) PR8GkAM t-ISING......... RESERENCS.ON *001 INITI.A IZATInN OF PARAMETFRS VFLAG=1 CVFLAG6O VECTOR VALUES INIIT= VECTOR VALUES PCH=l VECTOR VALUES PRSWI1 VECTOR VALUES EPSI=0.0001 VECTOR VALUES MXITER=100 VECTOR VALUES DELTA=1 VECTnR VALUES TSW=1 VECTOR VALUES RUN=1 VECTOR VALUES MAXJ=10 VECTOR VALUES CHKt-l MAIN0005 *001 MAIN006 *002 MAIN001 *003 MAIN0020 *004 MAIN0025 *005 MAIN0030 *006 MAIN0040 *OC7 MAIN0045 *008 MAIN0060 *009 MAIN0061 *010 MAIN0067 *011 MAIN0068 *012 0H I-J PROBLEM TITLE THE INITIAL RUN NUMBER MAY BE SPECIFIED.. IF DESIRED, BY PUNCHING THE STATEMENT.., THE INITIAL RUN NUMBER ISXXX," STARTING IN COLUMN ONE. BEGIN PR^NT FORMAT $1Hl,514,16HRQA RESULTS FOR,2C6,11H DISCIPLINE, 1 S20,,29HSFPT 1, 1965 VERSIONE_,RAQ-i / 1HQ*$,.NAME.NAME( 1 ) EXECUTE TODAY.l (DATE,YEAR) LOOK AT FORMAT t2C6+*$RUNP. UNP: _.._ _.. __ ___._ WHENEVER RUNP'E.$THE IN$.AND.RUNP( ).E.$ITIAL $tREAD FORMAT 1 $S25.T47*$,RUN WHENEVER RUN.GE.1003, RUN.0 __...PRNI_ EOR:MAT H, S1 4 4t4H.$R,_iUNJ_ _I,-_.S_. R DATEYEAR ZERO. ( TR NS...TRANS( 2390),,FIRST...F IRST( 2000),DEST.. DEST _ 1..2. 2 3 LAST.LAST12ZQL. __ __ _ ____ __ ___ _. _ MAIN0070 *013 _MA IN 071_ _ *_013 MAINO72 *314 MAIN0073 _._15 MAIN0074 *016 MATN0075 *016 MAIN0076. *017 MA IN 077 *018 MAIN0095 *019.MAIN 0096 *......... 9.......... __DAIA.. R.ED... A.AD. IN.A.A RATE...A.BLE. PREPARED..................... FXFCUTF D I STC. P MAT hI'1O 1 r * * * * * *.* * * *. * ** * *9 * * * * * *4 * * *, * * * * * * * * * * * * * * * * * * * * * * * * _, 9 9 9 9 9 9 + r, cw + 9 9 9

CALCULATION AND TEST OF THE NUMBER OF EN.TR-IBS IN THE _ TRANSITION TABLE. ERROR. IS A SYSTEM POST MORTEM DUMP RPIUT INE. THR)UGH MAIN3, FOR IR.,l 1.G.2000 _ IMIN3.WHENEVER TRANS I).NE.O., NTRANS = I WMHENE.VER CHK.E.-1".TRANSFER TO MAIN4 PRINT COMMENT $6 ****fTRANSITION TABLE EXCEEDS 2000 ENTRIES *1 $ -___________________ __ EXECUTE ERROR' CALCULATION AND TEST OF THE SCOPE OF THE PROB. VECTOR I___ _______________________________ A I N4 SCOPE=0..I IRInUH MATNS, FOR T=1,1,.G.NTNRANS SCOPET = DEST(I) + LAST(I). - FIRST(I) -_ _ ______ MWHENFFR LAS T(IT).G.CCnPFtSCOPF=I AST(I)__ MAIN5 WHENEVER SCOPET.G.SCOPE,.SCOPE = SCOPET WHENEVER SCrnPF IE.500F E TR:ANSFF R TO MATN6 PRINT COMMENT $6 ****THE LINEAR INDEX OF AT LEAST ONE STAT F.EXCEEDS 5000$ PRINT FORMAT $1HOS14,H+;ONE INDEX IS +.,I11*$.SCOPE ARATFM=I. ____ D iADD=. 'EXECUTF TnUTIr3 __ ___ CONVE'RGENT STOCHASTIC MATRI. PREPARED AND TESTED IF REQUIRED **~***********t******************************~*~******~ *****~ MAIN6 EXECUTE STOCAL. WHENEVER TSW.E.1, ExECUTTEJST.__. __ ____ CHOICE OF INITIAL VECTOR _ ENEVERI.E_____________________________ PRINT COMMENT $8 THE INITIAL ITERATE WAS DETERMI 1 NEn BY THE SUBROUTINE ESTIM.$ EXECUTE ESTIM. OR WHENEVER INIT.E.1. R- E.AU _EFORMAT IN IT1R, TEMP..TEMP( 6 )___ PRINT FORMAT INIT1,TEMP(1)....TEMP(6) DIMENSION TEMP(6) VECTOR VALUES INIT1=$.1H8,S14,38HTHE INITIAL ITERATE WAS THE R._ _1_iESULrL_ _ Ft_$.. _... _.. __..__ ___...__ ___ VECTOR VALUES INI TR=$S15,.6C6*$ MAINO111 MAIN0112 MAIN*112 MAIN0113 _MAIN114_ MAIN0115 *022 *023 *024 *024_ *025 01 01 001 01 MAIN0116 *026 MAIN0117 *027 MAIN0118 *028 MAIN*118 *029 MAINOl9 *030 MAIN0120 _ D031 MAIN0121 *032 MAIN0122 *032 MAIN0123 *033 MAINO?14 *034 MAIN0125 *035 MAIN0126 Q36 MAIN0140 *037 MAINO150 *038 MAIN0170__ * 039 MAIN0171 *040 MAIN0172 *040 MAIN0180 *041 MAINO090 _MAINO 191 MAIN0192 MAIN0193 MAIN0195.MANM0196 MAIN0197 * 042 *043 * 44 * 045 *046 *046 *047 01 01 01 oi" 01 01 01 01 01

EXECUTE BREAD. (SCOPEP) MAIN0202 *048 01 EXECUTE BREAD. V(YVFLAG 1 )...V(VFLAG,SCOPEP) MAIN0203 *049 01 OR WHENEVER -INIT.E.2 MAIN0210 *050 01 EXECUTEF ZERO. V(VFLAG.l)...V(VFLAG,SCOPE ),V(:CVFLAt1.)..*.V(.CVF MAIN0220 *051 01 1 LAG,SCOPE)) MAINO221 *051 vFLAG = 0 MAIN0222 *052 01 GCFLAG = 1 MAIN0223 *053 01 PRINT FORMAT $t1H8 S14,H+.THE INITIAL ITERATE IS GIVEN BELOW4+/1H-$q*$ MAIN0224 *054 01 READ AND PRINT DATA ' MAIN0230 *055 01 OTHERWI SE MAIN0231 *056 01 PRINT COMMENT $8 THE INITIAL ITERATE WAS THE RES MAIN0232 *057 01 1 ULT OF THE PREVIOUS RUNS MAIN0234 *057 END OF CONDITIONAL MAIN0240 *058 01 INITIAL VECTOR NORMALIZEDO_ ___________ __ ~Gi2 vcr7 VECSUM=O. THROUGH VEC1, FOR I=11,,I.G.SCOPE ' VECSUM=VECSUM+V(VFLAG, I ) THROUGH VEC2, FOR I=1,1I.G.SCOPE VVFL AG. I =V(VFL AC. /VCSUIIM M.AIN0260 MAIN0270 MAIN0280 MAIN0290 MAIN0300 *059 *060 * 061 *062 *063 01 01 ITERATION ROUTINE AND CONVERGENCE TESTING NITER=O MAIN0440 *064 MAIN1 MAXIX=CVGTST.(BEGIN) MAIN0450 *065 SAVEl=V(VFLAG,MAXIX)1} M A MAIN.460 *066, MAIN2 SAVE2=SAVE1 MAIN0470 *067 EXECUTE ITr-R. ( 1 ) SAVE =V (VFLAG, MAX IX) WHENEVER. ABS.LSAVE2-SAVE 1).LE.EPSI TRANSFER TO MAIN2 MAIN0480 MAIN0490... KI A _r- *068 *069 LVTr r,. 'q l - 1.Trn Tn Aq KC~ In t'n A" T I It * UK.iMAI hLK. LE.l IK I tKtI KAN-tRK IU MAINI[ M/A1NU)UU *UIU MAIN0510 *071 INTEGER INIT,MXITERII, ISCnPE, VFLAG, N ITER, MAXIX, MAXJ - MAIN0520 *072 1 CVFLAG.DELTAFIRST,.DEST,.LASTEXIT1,EXIT2,JCVGTST.,R~UN,OATE, MAIN0530 *072 2 YEAR,SCOPEP,NAME,TEMP,SCOPET,NTRANS,CHK,.PCH,TSW.PRSW,RUNPUl) MAIN0535 *072..._.EQ JUI VALE.NCE_.. (TRANS.) sNTRANS)TRANLST(.21 ),,CHK)....M - A _...__... _ MAN545 PROGRAM COMMON TRANS(2001),FIRST(2000),DEST(2000),LAST(2000), MAIN0550 *074._- PI......X.. HTER, EPSIl A, EL I SN A MAX I......... _ _ _______ MAIL0560..074 ERASABLE DUMMY( 100),VFLAG,,CVFLAG, SCOPE,NITEER.V((O...1)*5000), MAIN0580 *075 _1 ARATEM.VECSUM.EXITLEXIT2.EXIT3.I.J.DATEYEAR..DIADD MAIN0585 *075 END OF PROGRAM MAIN0600 *076

$COMPILE MAD,EXECUTE,DUMP,I/0 DUMP,UMIN'T 3BJECTPUNCH OBJECT STOC000O MA) (39 AUG 1965 VERSION) PROGRAM LISTIN... REFERENCES ON _ EXTERNAL FUNCTION _. ENTRY TO STOCAL. RATEM=0. _ THROUGH STOCi, F3, I=1,1,I.G.NTRANS WHENEVER RATEM.G.TRANS(I),RATEM=TRANS(I) STOC1 CONTINUE WHENEVER RATEM.GE.O. ARATEM =1./.99 DIADD=.31 TRANSFER TO STOC3 END OF CONDITIONAL ARATEM=(.ABS.RATEM)/.99 DIADD=1. STOC3 THROUGH STOC2, FOR I=1,1,I.G.NTRANS WHENEVER FIRST(I).NE.DEST(I) TRANS(I)=TRANS( I)/ARATEI OTHERWISE TRANS(I)=(TRANS(I)/ARATEM)+DIAD3 END OF CONDITIONAL STOC2 CONTINUE INTEGER I NTRANS,VFLAG,:VFLAG,DELTA,SCOPENI 'ER,FIRST, DEST, LA 1 STMXITERINITPCH,MAXJSCOPET EQUIVALENCE (TRANS,NTRAIS) PROGRAM CJMMON TRANS(2031),FIRST(2000),DEST(2000),LAST(2000)i 1 PCHPRSJ,MXITER,EPSIDE.TA,INIT,RUNTSW,NAME(1),MAXJ ERASABLE DUMMY(100),VFLAS,CVFLAGSCOPE NITERV({O...1)*5000), 1 ARATEM,VECSJM,EXIT1,EXIT2,EIT3,IJ,DATE,YEARDIADD FUNCTION RETURN END OF FJNCTION 0)0791 09/09/65 4 39 31*8 PM *001 srOCoo0010o *001 STOC0020 *002 SrOC0030 *003. STOCOo40 *004 srOCoo0050o *005 01 STOC0060 *006 01 SrOCOO. o o*007 STOCD0090 008 01 STOC0100 *009 01 STOC0105 *010 01 srocolO6 *011 01 STOC0110 *012 SOC0115. *013 srOC0120 *014 sroC0130 *015 01 s6OCOL46O 016 01 01 SrOCO15O *017 01 01 STOCO160 *018 01 01 SrOC0170 *019 01 01 STOC0180 *020 01 STOC0190 *021 STOC0200 ~021 Sr 0C0210 _ *022 srOC0220 *023 STOC0230 *023 sroC02-4 *024 STOC0246 *024 SrOC0260 *025 STOC0270 *026

$COMPILE MADEXECUTE,DUMP,I/O DUMP, IINT 3BJECTPUNCH OBJECT TESTODOO 0) 0791 09/00/65 4 39 37. 7 PM MAJ (9Y AUG L9b5 VRKSIUN) PROGRAM LISTINM......... ON 0\J REFERENCES ON EXTERNAL FUNCTION ENTRY TO TEST...******-'** ****~~~** **w**'**.'~~~~**,**-**-.** '"',**''.' TEST FOR LAST(If-FIST(I)=R*DELTA, WHERE R IS AN INTEGER ZERO.(V(CVFLAG,1)...V(CVFLAG,SCOPE)) EXITI1= THROUGH TESi, F3R I-=,',I-.G.NTRANS THROUGH TES2, FOR J=FIRST( I) tELTA,J.G.LAST(I) TES2. V(CVFLAG,J)=V(CVFLAG,J)+TRANS(I) __ _ WHENEVER J.E.LAST(I)+DE-TA.A,3.LAST(I.).LE.SCOPE.AND.FIRST(I) 1.G.O, TRANSFER TO TES1 EXITI=EXITl+l WHENEVER EXITI.E.1,EXECJTE TOUTT..EXECUTE TOUT1. TESI CONTINUE -.... TEST FOR STATE ROW SUMS - 1. AND POINT ROW SUMS = O. EXIT2=0 THROUGH TES3, FOR I=l,1,I.G.SCOPE WHENEVER VTCVFLAGi).E.D.. TRANSFER TO TES3 OR WHENEVER.ABS.(l-V(CVFLAG,I)).LE.l.E-7 TRANSFER TO TES3 OTHERWISE EXIT2=EXIT2+1 WHENEVER EXIT2.E.1,EXECJTE TDUTC. EXECUTE TJUT2. END OF CONDITIONAL TES3 CONTINUE TEST FOR POINT COLU4NS = 0 TEST0030 TEST0040 TEST0050 TEST0060 TEST0070 TEST0090 TEST0095 TEST0096 TESTO 100 TESTOlO TEST0120 *003 *004 *005 *006 *007 *008 *008 *009 *010 *OiO.011.*012 01 02 01 01 01 01 01 01 01 01 01 01 01 01 01 01 *001 TESTOO10 001 TEST0020 *002 TEST0130 *013 TESTO140 '014 TESTO150 *015 TEST0160 *016 TEST3170 *017 -TESTO180 *018' TEST0190 -o'9 TEST0200 *020 T=ST0210 *021 TEST0220 *022 TEST0230 *023 TEST0240 '-024.. _......... EXIT3=O THROUGH TES5, FOR I=1,, I.G.'TRANS THROUGH TES4, FOR J=DEST(I),DELTA,J.G.LAST(I)-FIRST(I)+DEST(I 1 ) TEST0251 TEST3252 TEST0253 TEST0254 *025 *026 *027 *027 01

WHENEVER V(CVFLAG,J).NE.0,TRAISFER TO TES4 EXIT3 = EXIT3+1 WHENEVER EXIT3,E. 1EXECUTE, TOUTM EXECUTE TOUT1. TRANSFER ro TES5 _________ TES4 CONTINUE TES5 CONTINUE WHENEVER EXIT1.NE.0.OR.EXIT2.SEO.OR.EXIT3. NE.0 EXECUTE 1 TOUT3. FUNCTION RETURN __ ________ INTEGER CVFLAG,EXIT1,I,'qTRAJS,JFIRST,DELTA,LASTEXIT2,SCOPE, 1 VFLAG,NITER,DEST,MXITER,INIT,PCH,MAXJEXIT3 EQUIVALENCE (TRANSNTRAMS) PROGRAM COMMON TRANS(2031),FIIST(2900),DEST(2000),LAST(2000), 1 PCH,PRS, MX ITER,EPSI,.DE TA, IlIT,RUN,TSW,NAME (),MAXJ ERASABLE DUMMY(100),VFLGS,CVFLAG^SCOPE,NITERV((O...1)*5000), 1 ARATEM, ECSUM,EXIT1,EXIT2,EXIT3, I J,DATE,YEARDIADD END OF FJNCTION TEST0255 *028 TEST0256 *029 TEST0257 *030 T-ST0258 *031 TEST0259 *032 TEST0260 *033 TEST0261 *034 T:ST0262 *035 TEST0263 *035 T=ST0264 *036 TEST0280 *037 TEST0290 *037 TEST0300 *038 TESTO310 *039 TEST0320 *039 TEST0340 *040 TEST0345 *040 TEST0360 *041 02 02 02 02 02 02 01 C7

ItrIMPTI F MAn. FXFrtlTF.nrIIMP I /nOnMP.PR:INT nRECT.PIUNCH OBJECT T[GUTUY0OO C0 Cs 2- Ij '9/2-n/65 1-0 43 37.5 PM MAO (09 AUG 1965 VERSION) PROGRAM LISTING......... ON R:EFERENCES ON_ E 'AL FUNCTION _........... EX T E_ N AL.__F_U__N............................................................._._._..... _....... ENTRY TO TOUTT. __ PRINT COMMENT $1. THE FOLLOWIN_G QUADRUPLE'S HAVE 1 AN ERROR IN FIRST OR LAST. J IS FIRST + N*OELTA, WHERE N IS 2 $_ _ _ _ _ _ _ _ PRINT COMMENT $ AN INTEGER.$ ___PRINT F___ MAT TI__ TLET _............... FUNCTION RETURN. ENTRY TO TOUTM. PRINT COMMENT $1 THE FOLLOWING QUADRUPLES HAVE 1 AN ERROR IN OEST. J IS A POINT INTO WHICH TRANSITION IS INDI _....-ATE D......... PRINT FORMAT TITLET ___..__._[-_ 0 T N...F ON RETURN. _..........-..-_ _ _...................._ _.._-_...........__. ENTRY TO- TOUT1. PRINT FORMAT OUT1, I,.TRANS( I),FIRST( I),OEST(I), LAST(I ),.J _........................................... C. R.. _ NTRY TO T U............TO..... _.T_.__._C~ PRINT COMMENT $1 THE FOLLOWING ROW SUMS, SUM (i) I. OF THE STOCHASTIC MATRIX ARE INCORRECT.$ PRINT FORMAT TITLFC...T N._-. TT. _ET. URN - -......................................... TENTRY O.. T OUT.2.................................._... PRINT FORMAT OUT2, I,.V(CVFLAG,I) FUNCTION RETURN._. -... ENIRY. TO._ TOUT3....... PRINT COMMENT $1 THE ABOVE IS A FATAL ERROR IN T 1 HE TRANSITION TABLE$ PRINT FORMAT $1HC,,SL4,H+IF ((U)) IS THE MATRIX DESCRIBED BY D _-__.L1.. ISC.EL.. _THENT_ I LE DSRIE-_E. THESCM__RLS+._ 4 ~ 7,H '_._( U) ) ',.E. 2 H'*((I))'*$,1./ARATFM,OIAnD PRINT FORMAT TITLED THROUGH TU, FOP, I=1,1,I.G.NTRANS __IU_............ PRINT OR.MAT OUT1, ITRANS.( I,.FLST.(_I).l,_DEST.I.)I._LAST( ) WHENEVER SCOPE G. 5),0 TR'ANSFER TO TUO PRINT COMMENT $1 THE FOLLOWING ARE THE ROW SUMS 1 OF THE STOCHASTIC MATRIX. PRINT COMMENT$8$ THROUGH. TUT, FOR. I=1i8,I+7.G.SCOPE. TUT PRINT FORMAT OUT4,I,.I+7,V(CVFLAG,I)...V(CVFLAGI+7)...........__...._!. C L.,._ _SCO..PE, PR I.E T FMA..O,_I!,_SCOPE4 V _CV, I_. )...V 1 (CVFL G,SCOPE) TUO EXECUTE ERROR. TOUTO01O ' 001 TOUT0020 * C02 TOUT0O30 *003 TOUTOC40 *003 TOUTOC'45 *003 TOUTOC;46 *. 004 TOUT0050 *005 TOUTOC60 *006 TOUT0061 *007 TOUT0O62 *008 TOUT0063 *008 TOUTO.A.64 *908 TOUTOC65 - 009.TOUTOp66.._*..._ 01.C TOUTO070 *011 TOUTO380 *012 TOUT0:9__ 0313 TOUTOIUO_ * 014 TOUT0110 *315 TOUT0120 - 015 TOUT313C *316 TOUTO.14..... 017 TnUTO.IS..... 08_ TOUT0160 *019 -..... T.U 1 70 * 020 TnUTCIPO 021 TOUTO190 *022 TOIUT 3195 022 TOUT02CC *'23 U T O02.1 _........._* 2 3_.... TOUT02C,2 *02? TOUT0220 '024 T rll T O 2 '~" * 02 4 TOUT 02?: *25 TnUTO24' - 026 TnUTC45S *027 TOUT2. * "'"8 TOUT C 26' *C28 TOUTOZ61 *029 Tn IJT 0 27 *03'. TOUTO2 8: *-':31 TOUT02~6 * 32 TDrJT 029- * C3: 01 01

VECTOR VALUES TITLET=$'H8.S15, 1HI,S8,8HTRANS (I).S 1I, 8HFIR9T( I 1 ) S8.7HOEST( S I S97HLAST(I),Sll tHJ/1H $ __..._. _ VECTOR VALUES OUT=II17,F15.8,4116*$ ___vE C TOR.VAL UES.. T.ITL EC=$ 1H 8 TS1..,. S9, UM.LI _.I._H __ __ _. VECTOR VALUES OUT2=$S14,lI4,S4,TFl1.8*$ _ VECTOR VALLES._I TLE-=$ 1HS,,S15_1-L. SH.T _-R.AA L.SL 11:,8HELRSLI 1 ) S8,7HDEST(I),S9,7HLAST(I)/1H *$ V ECTOR VA UES nUlL4-i S14,I4.3H -. T 4,S3. 8 ( F8. S3 ) / H *$ ____NlTEGR. F E. RSTIIDI1ESTL, LASTIT_,IJ_ VF A Z NZTRANSVFELAGn ELTA.S COPE.E_..1 NITER,MXITER, INIT,PCH,MAXJ,EXIT1,EXIT2,EXIT3....E..QU 1.V ALEMCE- _lR ANSsN.. I RANS. _............-.. PROGRAM COMMON TRANS(2001 ),FIR:ST(2000),DEST.(2000),LAST(2900), 1 PCH ~ PRSW,.MX I TER?FP S. I tEL TA, I N I T, RUN, TSW,.NAME ( 1 ).,MAXJ ERASABLE DUMM ( O) VFLAGtCVFLAG SCOPE,NITER ( ( 0... 1)*5300 ) 1_ARATI.EM. V~CSUI X.E XIIL F X IT2 E-,.XXI3 I 3-.IL.DAIE_4YF-AR L _ DID__A __DD_. FUNCTION RETURN..... E D_..E._.FUNCTI N... -.-..... TO.UT0380 *034 TOUT039.' *034 TOUT0409 *035 TOUT0410: *036 TOUT0420 *037 ___IIOUI.T43..* - *33 8. TOUT0440 *038 TOUT046O *039 _..ILQUIQ49 _.*.... 4.:C -4 TOUT0500 *040 T... T 5IO..............._41. TOUT0520 *042 _ TOUT0530 *042 TOUT0540 *043 __ _TUT_0.5-... —__-._ -34 3 _ TOUT0590 *044 _ T.OUTIQ6.00-..._5 N0 0O

$COMPILE MAD,EXECUTE,DUMP, I/O DUMPPRINT OBJECT,PUNCH OBJECT ESTMOOO0 003840 02/04/66 2 15 17.4 PM MAD (03 JAN 1966 VERSION) PROGRAM LISTING........ REFERENCES ON *.C0 EXTERNAL FUNCTION ESTMOO1O *001 ENTRY TO ESTIM. ESTMO020 *CC2 LAMSTR=O. ESTM0030 *CC3 MUSTR=O. ESTM0040 *C04 T2=0 ESTM050 * 005 T3=0 ESTM0060 *006 THROUGH ES1, FOR I=l,l,I.G.NTRANS ESTM0070 *007 DIFSTR=FIRST(I)-DEST(I) ESTM0080 *008 Tl=((LAST(I)-FIRST(I))/DELTA)+1 ESTMOO90 *09 WHENEVER DIFSTR.E.O,TRANSFER TO ES1 ESTMO100 *010 WHENEVER DIFSTR.G.O,MUSTR=MUSTR+TRANS(I)*T1 ESTMO.11Q *011 WHENEVER DIFSTR.L.,LAMSTR=LSTR STR+TRANS(I)*T1 ESTM0120 *012 T2=T2+.ABS.DIFSTR*T1 ESTMO130 *013 T3=T3+T1 ESTMC140 *014 ES1 CONTINUE ESTMO150 0C15 RHOSTR=(LAMSTR/MUSTR).P.((l.*T3)/T2) ESTM0160 *016 'm~**~* **1~" " 4F*********m~i * 9\O^~~ ~ESTIMATE ONLY VALUES WHICH ARE LARGER THAN DESIREO ACCURACY WHEN NORMALIZED. Tl=l ESTMO170 *017 T2=SCOPE ESTMO180 *018 ESCOPE=.ABS.(ELOG.(EPSI)/ELOG.(RHOSTR)) ESTMO190 *019 WHENEVER ESCOPE.GE.SCOPE,TRANSFER TO ES2 ESTM0200 *C20 ~WHENEVER RHOSTR.GE.1., T=SCOPE-ESCOPE ESTM021Q *021 WHENEVER RH0STR.L.1.,T2=ESCOPE ESTM0220 *022 EXECUTE ZERO.V(1)....V(SCOPE)) ESTM0230 *0C23 ES2 RHON=l. ESTM0240 *024 CVFLAG=1 ESTM0250 *025 VFLAG=O ESTM0260 *026 THROUGH ES3, FOR I=T1,1,I.G.T2 ESTM2.70 *027 RHON=RHON*RHOSTR ESTMO280 *028 ES3 V(I)=RHON ESTMO290 *029 ZERO OUT VALUES AT POINTS BY PERFORMING AN ITERATICN EXECUTE ITER.(1) ESTMO300 *030 FUNCTION RETURN ESTM0310 *C31 INTEGER T1,T2,T3,ESCOPE,DIFSTR,I,NTRANS,FIRST,DEST,LAST, ESTM0350 *032 1 DELTA,K,J,CVFLAG,VFLAG,SCOPE,NITER,MXITER,INITPCH,IAXJ ESTMC360 *032 EQUIVALENCE (TRANS,NTRANS) ESTMC370 *C33 01 01 01 01 01 01 01 01 01 01

PROGRAM COMMON TRANS(2CO1),FIRST(2 EST(T(2CCO),LAST(2000), ESTM0380 *034 PCH,PRSW,MXITER,EPSI,DELT A, IN T,RUN,TSW,NIAME(1),MAXJ ESTM0390 *034 ERASABLE DUMMY(100),VFLAG,CVFLAG,SCOPEtNITERV((0...1)*5000), ESTM0410 *035 ARATEM,VECSUM,EXIT1,EXIT2,EXIT3,.,J,OATE,YEAR,CIACC ESTMO420 *C35 END CF FUNCTION ESTM0430 *036 0

SASSEMBLEEXECUTE,DUMPPRINT OBJECTPUNCH OBJECT ITEROOOO O01196 09/13/65 12 58 25.8 PM PAGE 1 (MO0 05 MAY 1965) UNIVERSITY OF MICHIGAN 13 SEP 1965 -DRnFriAM I:UTU InrTA I f: rr.T n I rT nucT:ACAaDl: ll::rt ARlI nir- f-i h — 1 --- nT fTFRE tlfl NIuTIAL TRANSFER VECTOR 00000 INITIAL.TRANSFER VECTOR. 00000 TRANSFER VECTOR LENGTH TRANSFER VECTOR LENGTH' LUWC41 CKA4ADLE UEr&NCU 144UO Iu a *I %wuw PROGRAM ENTRY POINT LOCATIONS AND NAMES. 00000 I TER PROGRAM LISTING (NUMBERS IN OCTAL) H ENTRY ITER 00000 0634 00 4 00060 010 IT i- TE- X --- XR4,4 00001 0500 60 4 00001 GO CLA* 1,4 00002 0601 00 0 00064 010 STO N 00003 0500 00 0 27511 CO CLA DELTA OOO"4 0767 00 0 00022 00 -ALS 18 00 05 0622 00 0 00044 010 STD D1 00006 0622 00 0 00045 01 —0 STD -— 0- 2 * 00007 0534 00 1 77632 010 START LXA VFLAG,1 00010 0500 00 1 00063 010 CLA AD2,1 00011 0621 00 0 00042 010 STA RV1 00012 0621 00 0 00043 010 STA RV2 0013- 0621 00 0 00025-.O STA. RV3 - - 00014 05000 0 0 77632 010 CLA VFLAG -0015 0322 00 0 00065 010 ERA =1t 00016 0734 00 1 00000 CO PAX,1 00017-500 00 1 00063 C10 A TO CLAT --- AD2-' 00020 0621 00 0 00040 010 STA OPV1 00021 05000 00 77627 010 - CL-A -- - N-I-TER-. 00022 0400 00 0 00065 010 ADD =1 00023 060 0 0 77627 010 STO NITER 00024 0534.00 1 77630 010 LXA SCOPE,1 00025 060 00 1 000000 00 RV3 STZ **,1 00026 2 00001 1 00025 010 TIX - 1,1,1 <* 00027 0534 00 -4 13721 00 LXA TRANS,4 00030 0500 00 4 27504 00 KI CLA -LAST-4 00031 0767 00 0 00022 00 ALS 18 00032 —0622- 0 0 00037 010 -STD- - - - 00033 0500 00 4 17642 00 CLA FIRST,4 -00034 0734 0 00 '0 2 0 O. PAX y - 2 00035 0500 00 4 23563 00 CLA DEST,4 00036 0734 00 1 30000.00 PAX --- - 00037 3 0000- 00A-O-60 —0 - l T-XH L, 2,* 00040 0560 00 2 00000 00 OPV1 LDQ **,2 00041- 0260 -0 4-r372 -72 - -' F — M- P - ~TRANS,'4 00042 0300 00 1 00000 00 RV1 FAD *,1 00043 0601 00 1 00000 00 RV2 STO *-, 00044 1 00000 1 00045 010 D1 TXI *+1,1,** 00045 I 00000 2 00037 022 TXI — L-,2,**00046 2 00001 4 00030 010 L2 - -- TI-X K 1,4, 00047 0500 00 0 77632 010 CLA - VFLAG (ninf n'Ar l ri. n 77At1 /n)i CTn rU- Ar DO ITERO01 ID ITER0002 VVvjV Voul vi v (FOOL 00051 0322 00 0 00065 U 1U 010 a IU t, VrL A — CRA - - = 1- —.-. ID ITER0003

PAGE 2 00052 0601 30 0 77632 010 g L —'US~-' D5 — O0 6-4- G IU'00054 0402 00 0 00065 010 00055 0601 00 0 00064 0 10 -00056 0520 0 0 00064 010 '0057 -"0020 CO 0 00007 010 00060 0774 00 4 00000 00 -00061 2 O 0Z o 4.. oo..-.. 00062 0 03.00 0 77626 010 UNIVERSITY OF MICHIGAN 13 SEP 1965 STO VFLAG............ CL...........N............. SUB =1 STO N ZET N TRA START XR4 AXT **,4 TRA 2,t4 AD1 PZE V 00063 O — O 000 0 066016 010 AD2 PZE V-5C00 00064 0 00000 0 00000 00 N PZE 13721 TRANS PGMCOM 2001 17642 FIRSI PGMCOM 2CCO 2356'3 — -— " -DEST -PGMCOM 2C000 27504 LAST PGMCOM 2CO0 27510 PGMCOM 3 PCLIST DELTA 27511 DELTA 27517 PGMCOM 5.................... -.... 77632 ERAS 100 ERL1ST VFLAG,CVFLAG,SCOPE,NITER 77632 VFLAG 77631 CVFLAG 77630 SCOPE 77627 NITER 77626 V ERAS 100CO END PROGRAM LIT~RALS 00065 OOOOGOC301 CO **** 00066 IS THE FIRST LOCATION NOT USED BY THIS PROGRAM. REFERENCES TO SYMBOLS DEFINED IN PROGRAM. TYPE SYMBOL VALUE REFERENCES TO SYMBOL (SYMBOL TYPES - R - RELOCATABLE, A- ABSOLUTE, E - ERASABLE, T - TRANSFER VECTOR) R AC1 0006a 30062 -R A02.-.. O063,3T01, 00 17, C0063 E CVFLAG 77631 00050,0C0065 R 01 00044 00C05,03'J44 R D2 00045 006,00()45 A DELTA 27511 00003,00065 A DEST 23563 0CG35,00065 R FI' 0: 041 "00'' 41 A FIRST 17642 C0033,00065 R ITER 00000 OC!' '^COC, 0 0 R KI C0n030 00030,00046 R LI 0' ' 3 7 Cr'. 32,0; 037, 00045 R L2 00046.0037,00046 A LAST 27504 n'O 330, 0^65 R N 0' 64 "0 I2,'r, 53, -0055,00i356 00064 E NITER 77627 CO0: 21,.' ~C;23, G 65 R OPVl 0004C C 4 20, 2 0r:'43 R RV1 00042 )00011, 90 042 R RV2 00043 30012,0 43 R RV3 G ' 0 25 C7) 13, O C25 SCOPE 7763C C0024,00,'65 R START,':? 7 00137, C057

PAGE 3 A TRANS 13721 0CO27, '30041, &C065 E " V - 77626 S00 ^,-2',C6. ',..O65 E VFLAG 77632 0'007, O, 014, 0"O47, O' 052,00065 R XR4 006 6 tOCO, C0,O;06Q REFERENCES TO PROGRAM LITERALS LOC REFERENCES 000 6 '0-O5,'00022,000 51, OCO 54 FLAGGED LINES FATAL 0 NON-FATAL 0 MACRO CALLS CORRECT 0 INCORRECT 0 UNIVERSITY OF MICIIGAN 13 SEP 1965 **** DEFINITICNS ***.* CARD COUNTS SYMBOL 26 INPUT 71 CPERATION 0 CUTPUT 4 MACRO 0 RMT 0

ArnMPTI:F MADn. FYtIITE..nMP' T/n n DMP.PRINT O EnRIF-T.PINH ORJECT0 002300 09/20/65 10 43 46.4 PM MAD (09 AUG 1965 VERSION) PROGRAM LISTING......... _EXTFRNAL FUNCTtI N ENTRY TO CVGTST. IF NITER - O0 SET TRAP PROCEDURE SO THAT UNDERFLOWS _..._ PRnODUCE ZERnS., ALSIn ElN) _liSTTARIF iN!rFX FOR: WFAK TESaT___.... _. ***~~~*~~~~~~~~~e*~~~~~ * *~***~~~;~~*~~*;~******~**.__ VFFTnR VALUES FTSW=, WHENEVER FTSW E.O WHFNFVFR NTTFR'.GE, MXIT_ __TR _ ____ MAXDIF=O. _TRANSFER TO CV3 _____ END OF CONDITIONAL FTSW=1 EXECUTE FTRAP. T.IlROU L=J C VO. V ( VFA,. DOR. I.F. SCOPE __ FUNCTION RETURN I _. —ED..E - O..._ IT OND N T_ ONAL ---_-_ ____ *001 — C VGT*n001 ~01 CVGT0002 *002 CVGT30Q3 *OC3 CVGT0004 *004 _V G.TIO05..__._.*05. CVGTOO06 *006 CVGTOOQ7 *_.___.7 CVGT0008 *008 CVGTOO14 *009 CVGTOO15 *010 CVGTOO16 _ *011 CVGTOO17 *012 CrTQO 18 __*13. 01 02 02 02 01 01 01 01 --- FTNn MAXIMIIM DIFFFDfFNCEf AN TCFT ArcATNST FPST MAx IF= O. THROUGH CV2, FOR I=1,1, I.G.SCOPE _. _fD I _ABS..(V VFLA_ ) - V( CF L A GI )) WHENEVER MAXDIF.G.DIF, TRANSFER TO CV2 __ M AXDiF = IF... _... __ MIX=I TONT INUIF.....E..EVER..:E.E._MAXDIF_....LE_. E.PS _PRINT COMMENT $8 ***-**************j**-*********** SUCC _ T L.OESSL_ LUN *****************************$ _ __ C_ TRANSFER TO CV1 CVGT0O30 *014 CVGTOO40 *015 CVGTQ050Q __ *016 CVGTO060 *017 CVGT0O070 _ *018 CVGT0080 *019 CVGTOC90 *020 CVGTO100 *021 CVGTO105 *022 CVGTO110 *022 CVGTO115 *023 CVGTO120 *024 CVGTO125 0__ *25,. CVGTO130 *025 01 01 _.0.Q1..... 01 01 01 01 01 01 OR WHENEVER MITER.GE. MXITER CV3._.... PRI.I.CMME.N_ $ 8**.**_**. *__THE SPECIFIED NUMBER OF ITERATIONS_ 1 HAS BEEN SATISFIED ********** $ ESTIMATE SUBDOMINANT EIGENVALUE AND ERROR.

CMl LDIF=MAXDIF.P. 1./NITER-) WHENEVER LDIF.GE.1..OR.LDIF.L..5 PRINT FORMAT CVFL,,EPSI, MXI TER,NITER,MAXDIF MAXDIF O. OR WHENEVER. 1.-LnlF).P. (1./NITER).L.LDIF ____ PRINT FOR:MAT CVFL, EPSI:,MXITERNITERMAXDIF O., LDIF OTHERWISE__ __ PRINT FORMAT CVFL,tEPSI,.MXITER,.NITERMAXDIFtMAXDIF/LDIF,LDIF END OF CONDITIONAL RESET TRAP PROCEDURE SO THAT UNDERFLOWS AGAIN CAUSE ERROR: REiTRNS. ~ ***~~*****~ ****** #~~~~ *; *~>~*** ~ ~ * EXECUTE NTRAPw FTSW=O _EXECUTE OUT 2:. CVGT0131 *026 01 CVGTO132 *027 01 CVGT0133 *028 02 CVGT0134 *029 02 CVGT0135 *030 02 CVGT0136 031 02_ CVGT0137 *032 02 CVGT0138 *033 02 CVGT0139 CVGT0140 CVGTO 145 *034 *035 *036 0.1 01 01 THE ER:ROR RETURN TR:ANSFERS TO 'BEGIN' IN THE MAIN PROGRAM ERROR RETURN_____.___ OTHERW I SS _____ _ ______ FUNCTION RETURN MIX END OF CONDITIONAL VECTOR VALUES CVL=$1H-_5S1435_HCONVERGFNCE TBEST CRITERION__EP 1 SI) =F11.8/lHO,,S12,37HMAXIMUM ALLOWED ITERATIONS (MXITER) =15 2 /1HO.S31.,18HITER'ALONS_._T_ AKEN = Il5LHO, S3.4HMAXI MUM _..DFFENC. 3 E BETWEEN LAST TWO ITERATES =E8.1/1HO, S23,26HESTIMATED ORDER 4 OF ERROR =E8.1.S3.33H(APPROXIMATF CONVERGENCE FACTOR =F9,7, 5 1H)*$ INTEGER I,SCOPE,VFLAG,CVFLAG,MIX,DELTANITERtFIRST-OEST,LAST, lMX LTERtIl. _ _ ___ -- - PROGRAM COMMON TRANS(2001),FIRST(2000),DEST(2000),LAST ( 2C.OC), 1 PCH,PRSW.MXITEREPSI,DELTA,.INITRUN.TSW,NAME(l1),MAXJ - _ ERASABLE DUMMY(1OO),VFLAG,.CVFLAG,SCOPE,NITERV ( (...1)*500C),. 1 ARAIEMVEC. SU E.XI_ 1,EXL2 F E X_ 3TJ3_SL JDTE ARaDI ADD. END OF FUNCTION CVGTO150 *03._._ 01 CVGTO155 *038 01 CVGT0160 *039 01 CVGT0165 *040 01 CVGTO170 *041 CVGT0180 *041 CVGT018 O2__._ *041 CVGT0184 *041 CVGT0186 *041 CVGTO180 *041 CVGT0185 *042 CVGT0186_._ *Q42. CVGTO190 *043 CVGTO200 * 943 CVGT0205 *044 CVGTO21 ~ *044 CVGTG90 *045

$COMPILE MADOEXECUTEDUMP,I/O DUMP,PRINT OBJECT,PUNCH OBJECT' OUTTOOOO MAD (09 AUG 1965 VERSION) PROGRAM LISTING......... REFERENCES ON EXTERNAL FUNCTION ENTRY TO OUT2. WHENEVER PCH.NE.1,TRANSFER TO UT9 PRINT COMMENT S- THE SOLUTION VECTOR HAS BEEN PUNCHED$ OUTPUT VECTOR CARDS WILL HAVE THE FOLLOWING ID IN COLUMNS 73-80 (1) MODEL NAME 12) DATE (3) RUN NO. XX 001 ALL SUCCEEDING CARDS WILL HAVE RUN NO. XX VECTOR CARD NO. * *.**.* *** **.** *.*** *.* *.**.** I.*.****** *** *************** * ***m * ** PUNCH FORMAT S1H,S14,2C6,S1I4HRUN,13,2H,,2C6,S23,C6,C2*S, I NAMENAME(1),RUN,DATEYEARNAME,NAME(1) EXECUTE NSDBCD.(DATES S) EXECUTE BPUNCH. SCOPE) EXECUTE NSDBCD.((BNBCD.(RUN).LS.18).V.$OO 0$,$01s) EXECUTE BPUNCH.(V(VFLAG,l)...V(VFLAfG SCOPE)I UT9 WHENEVER PRSW.E.1,EXECUTE PRQUT. RUN = RUN + 1 INTEGER,NITER,MXITER,I,SCOPE,J,K,DELTAVFLAGCVFLAG,FIRST,DES 1 TLASTINIT,MAXJ,S,RUN,DATE,YEARNAMEBNBCD.,PCHPRSW PROGRAM COMMON TRANS(200/1),FIRST(2000),DEST(2000),LAST(2000), 1 PCH,PRSW,MXITEREPSI,DELTA,INIT,RUN,TSW,NAME'(1) MAXJ ERASABLE DUMMY(100),VFLAGCVFLAG,SCOPE,NITER,V((0...1)*5000), 1 ARATEM,VECSUMEXITl,EXITZEXIT3,I,JDATEYEARDIADD FUNCTION RETURN END OF FUNCTION 001632 09/15/65 6 18 16.2 PM *001 OJTTOO10 001 _ OJTT0020 *002 OJTT0300 *003 OJTT0320 *004 — 4 G\ OJTT0351 OJTT0352 OUTT0354 OJTT0355 OJTT0356 OUTT0357 OJTT0370 OJTT0380 *005 *005 *006 *007 *008 *009 *010 *011 OUTT0430 *012 OJTT0440T * 012 OJTT0470 *013 OJTT0480 ~*313 OJTT0495 *014 OJTT0496 *014 OJTT0510 *015 OUTT0520. '*01.6

SCOMPILE MAD,PRINT OBJECT,PUNCH OBJECT POUTOOOO 002377 09/21/65 3 42 17.3 PM MAD (09 AUG 1965 VERSION) PROGRAM LISTING........ REFERENCES ON O001 EXTERNAL FUNCTION PJUTOO10 O 001 ENTRY TO PROUT. P3Ur0020 *002 PROBABILITY VECTOR PRINTED OUT PRINT FORMAT $1H1,S20,H+THE FOLLOWING ARE THE LIMITING STATE PDUT0060 *003 1 PROBABILITIES, PROB (IJ).+/S20, H+ PROB(I) IS THE SUM OF PR POUTO070 *003 2 OB(I,J) OVER ALL STATES WITH THE VALUE I FOR THE FIRST VARIAB POUT0080 *003 3 LE+/IH-*$ POUT0090 *003 WHENEVER (SCCPE/MAXJ)+1.G.2500.OR.MAXJ.G.25CO PJUT0095 *004 PRINT COMMENT $ *** MAXJ, THE MAXIMUM COLUMN INDEX HAS B POUT0100 *005 01 1 EEN SET TO FIFTY TO ALLOW MORE EFFICIENT PRINTING CF THE PROB POUT0105 *005 2 ABILITY VECTOR.$ POUT0106 *005 MAXJ=50 POUTOIlO *006 01 END OF CONDITIONAL POUT0115 *007 01 PRINT FORMAT $1H-,S3,1HI,S4,8H PROB(I),S40,9HPROB(I,J)/1H-, POUT0120 *008 1 S19,3HJ =,S2,10(11,S1O)*$( =O,1,I.G.9.OR.I.G.(MAXJ-1),I) PDUT0121 *008 VECSUM=O. POUT0122 *009 EXECUTE ZERO.(MARG1...MARGI((SCOPE/MAXJ)+1),MARG2...MARG2(MAX POUT0127 *010 1 J)) POUT0128 *010 THROUGH UTTA, FOR K=1,MAXJ,K.G.SCOPE POUT0130 *011 THROUGH UTTB, FOR S=K,1,S.G.(K+MAXJ-1).OR.S.G.SCOPE P3UT0140 *012 I=(S-1)/MAXJ POUT0150 *013 dJ=S-1-I*MAXJ POUT0160 *014 MARGl(I)=MARG1(I)+V(VFLAG,S) P3UT0170 *015 UTTB MARG2(J)=MARG2(J)+V(VFLAG,S) POUTO180 *016 VECSUM=VECSUM+MARG1(I) POUT0190 *017 WHENEVER (K+MAXJ-1).G.SCOPE POUT0200 *018 PRINT FORMAT LINE,I,MARGl(I),V(VFLAG,K)...V(VFLAG,SCOPE) PJUT0210 *019 01 UTHERWIS[ POUT0220 *020 01 PRINT FORMAT LINE,I,MARGI(I),V(VFLAG,K)...V(VFLAG,K+MAXJ-1) POUT0230 *021 01 U TA END GF CONDITIONAL POUT0240 *022 01 VECTOR VALUES LINE=$1HO,14,S2,FO1.8,S2,1C(F10.8,S1)/( P3UT0250 *023 1 S19,10(F10.8,S1))*$ POUT0260 *023 PRINT FORMAT $lH8,S15,H+ PROB(J) IS THE SUM OF PROB(I,J) OVER POUr0270 *024 1 ALL STALES WITH THE VALUE J FOR THE SECOND VARIABLE+/1H-, POUT0280 *024 2 S8,3HCHFCKSUM,S40,7HPROB(J)/iHO,S19,3HJ =,S2,10( 1,S10)*$, PJUTO290 *024 3 (1=0,1,I.G.9.OR.I.G.(MAXJ-1),I) POUT0291 *024 PRINT FORMAT $1HO,S6,F10.8,S2,10(FIC.8,Sl)/(S19,10(FO1.8, POUT0310 *025 1 Sl))*I,VECSUM,MARG2(0)...MARG2(MAXJ-1) POUT0320 *025 INTEGER NITER,MXITER,MAXJ,SCOPE,K,I,J,VFLAG,CVFLAG,S POUT0330 *026 DIMENSION MARGI1(2500),MARG2(2500) POUT0340 *027 EQUIVALENCE (MARG1,FIRST),(MARGI(2001),TTRANS),(MARG2,LAST) POUT0345 *028 1,(MARG2(2C01),CEST) P3UT0346 *028 PROGRAM COMMON TRANS(2001),FIRST(2CCO),DEST(2000),LAST(2000), POUT0350 *029 01 02 02 02 02 01 01 01 01 01 01

1 PCH,PRSW,MXITER,EPSI,DELTA,INIT,RUN,TSW,NAME(,MAXJ POUT0360 *029 ERASABLE DUMMY(100),VFLAG,CVFLAG,SCOPE,NITER,V((O...1)*5000), POUT0380 *030 1 ARATEM,VECSUM,EXITl,EXIT2,EXIT3,I,J,DATE,YEAR,OIAOD P3UT0390 *030 FUNCTION RETURN POUT0400 '031 END OF FUNCTION POUT0410 032

SCOMPILE MAD,PRINT OBJECT,PUNCH OBJECT ITEROCOO OOOD787 09/09/65 2 09 0.9 PM MAC (09 AUG 1965 VERSION) PROGRAM LISTING...... REFERE NCES ON ITERu010 *001 EXTERNAL FUNCTION (N) ITER0220 *001 ENTRY TO ITER. ITEROj30 *002 NORMAL MODE IS INTEGER ITER0040 *003 FLOATING POINT V,TRANS ITER0,50 *004 THROUGH END, FOR CNTR=1,1,CNTR.G.N ITER0060 *005 EXECUTE ZERO.(V(CVFLAG,1)...V(CVFLAGSCOPE)) ITER0070 *006 NITER = NITER+1 ITERO80 *007 THROUGH KLOOP, FOR K=1,1,K.G.NTRANS ITER0090 *008 J=FIRST (K) ITER0100 *C09 I=DEST(K) ITERO110 *010 CC]MPLP V( CVFLAG, I ) =V( VFLAG, J )*TRANS (K)+V(CVFLAG, I )ITER0123 *011 I=I+OELTA ITERO130 *012 J=J+DELTA ITEROi40 *013 KLOOP WHENEVER J.LE.LAST(K),TRANSFER TO COMPLP ITERO150 *014 VFLAG=VFLAG.EV. 1 ITERO160 *015 END CVFLAG=CVFLAG.EV.1 ITER0170 *016 PROGRAM COMMON TRANS(2001),FIRST(20'0O),DEST(200J),LAST(2000), ITERO130 *017 1 PCH,PRSW,MXITER,EpSI,DELTA,,INIT,RUN,TSW,NAME(1),MAXJ ITERO190 *017. _. ERASABLE OUMMY(!Cl),VFLAG,CVFLAG,SCnPE,NITER,V( (0... 1)*5000), ITER0210 *018 \0 1 ARATEM,VECSUM,EXITi,EXIT2,EXIT3,I J,DATE,YEAR,DIADD ITER0220 *018 EQUIVALENCE (NTRANS,TRANS) ITER0230 *019 FUNCTION RETURN ITER0240 *020 END OF FUNCTION ITER0250 *021 01 01 01 02 02 02 02 02 02 01 01

APPENDIX B SOME SUBROUTINES FOR USE IN GENERATING QUADRUPLES IN RQA-1 James H. Jackson B.1 PURPOSE A three-entry subroutine has been written in an attempt to simplify the programming of discipline subroutines for use with RQA-1. The entries and their functions are as follows: SETUP - Initializes parameters internal to the subroutine QUAD - Produces the quadruples which represent transitions between states (the off-diagonal elements of the Q matrix) DIAG - Produces the quadruples which represent the diagonal elements of the Q matrix. B.2 CALLING SEQUENCES (MAD) SETUP. (STATE.) STATE is the name of an integer-valued internal function whose value is the state number represented by the current values of the state variables. Its values may range from 1 through 5,000. QUAD. (RATE, CHANGE, REPEAT) RATE is a floating point expression whose value is the transition rate for the generated quadruple. CHANGE is an integer expression whose value is the amount by which the state number is changed when a transition indicated by the generated quadruple occurs. REPEAT is an integer expression chosen such that the generated quadruple will apply to STATE., STATE. +DELTA,..., STATE.+REPEAT*DELTA. DIAG. This entry has no arguments. 80

B35 RESTRICTIONS SETUP must be called before either QUAD or DIAG is called, otherwise an error comment will be printed and control will be returned to the system. MAXJ must be initialized before SETUP is called. QUAD must be called as many times as needed to produce all of the quadruples which represent transitions between states. The last call to QUAD must occur before DIAG is calledo After DIAG is called, no more calls should be made to the subroutine until RQA-1 has processed the generated data. Bo4 DISCIPLINE PROGRAM STRUCTURE The discipline program which calls SETUP, QUAD, and DIAG may have the following structure: EXTERNAL FUNCTION ENTRY TO DISCPLo READ AND PRINT DATA MAXJ = DELTA 3 -- SETUP. ( STATE ) THROUGH S. FOR --- THROUGH S, FOR -- WHENEVER -—, TRANSFER TO S WHENEVER --—, QUAD.( ----, ---, -— ) WHENEVER -—, QUAD. ( —,-,..) DIAGo FUNCTION RETURN FLOATING POINT EPSI, - NORMAL MODE IS INTEGER PROGRAM COMMON DUMMY (8004), PCH, PRSW, MXITER, EPSI, 1DELTA, INIT, RUN3, TSW, NAME(1), MAXJ VECTOR VALUES NAME ~= ---- INTERNAL FUNCTION STATE. (X) = END OF FUNCTION When the discipline subroutine is called, it reads a set of data. This data should include values of the parameters which describe the model as well as the values of PCH, PRSW, MXITER, EPSI, INIT, RUN, and TSW which are to be different from the preset values determined by RQA-k!.

After the data is read, MAXJ and DELTA are initialized. Then SETUP is called to initialize internal parameters to the subroutine. A set of nested interations is then begun. Each iteration variable is a state variable. Hence, within the innermost loop, all possible combinations of state variables will be produced. Since impossible combinations of state variables may not always be eliminated by constraining the set of values each variable may assume, the first conditional statement shown specifies a transfer to the end of the loop when any constraint of the model is violated. The remaining conditional statements generate quadruples which represent transitions between states whenever transstates whenever ransons may occurfrom the sta represented by the current values of the state variables. When all of these quadruples have been generated, the diagonal elements of the Q matrix are generated, and control is returned to RQA-1. B 5 EXAMPLE As an illustration of the programming of a discipline routine using SETUP, QUAD, and DIAG, one might consider a second-order Erlang channel with a finite queue: --- ( —.~- 2 -2 — N = Maximum length of queue = Arrival rate p = Processing rate He may choose to define the state variables as follows: A = Number of Erlang channels busy B = Number of Erlang channels in seco nd phase C = Queue Length The ranges of these state variables are the following: O<A < 1 O<B <A O<C <N Since RQA-1 will print a two-dimensional probability array, the state space must be described in terms of two state variables to facilitate interpretation of the results. One may choose one of these two variables to be the state of

the service and the other to be the queue length. Since the state of the service is completely defined by A and B, it will be denoted by S(A,B), whereas the queue length is simply C. S(A,B) may be defined by the following table: A B 0 1 1 0 0 1 S(A,B) 0 1 2 from this table, one may derive the following algebraic expression for S(A,B): S(A,B) = A + B Since S(A,B) has three the function STATE must be be expressed as follows: possible values, and since the minimum values of at least one, the value of STATE, S(A,B,C), may S(A,B,C) = S(A,B) + 3C + 1, or S(A,B,C) = A + B + 3C + 1. We may now define all possible transitions for the model by the following table: Condition O < C < N A = 1 B = 0 Event B+ Change 1 Rate 2|i I 1 < C < N C + -4 2[ B = 1 B+ C = 0 A -2 2t B = 1 BJ C 0 A + 1 X A =0 O < C <N-l C+ 3 X A = 1 85

From this table, one can infer that many of the transitions are repeated for the various values of the queue.length C. Consequently, C should be chosen to be the first index for the probability array, and DELTA should be chosen to be the number of values of S(A,B), or 3. MAXJ is then 3, By applying the above table of possible transitions and the program form of Section 4, one may obtain the following discipline program for the model: EXTERNAL FUNCTION ENTRY TO DISCPL. READ AND PRINT DATY MAXJ = 3 DELTA =3 SETUP.(STATE.) THROUGH S, FOR A = THROUGH S, FOR B = THROUGH S, FOR C = WHENEVER C. G. 0. t WHENEVER C. E. 0. i WHENEVER C. E. 1.o WHENEVER C. E. 0.o WHENEVER C. E. Oo 1 WHENEVER C. E. 0.O CONTINUE DIAG. A 0, 1, 1 0,1,] 0, 1,( ANDo A. AND o A. AND. B. AND. B. AND. A. AND. A A. G. 3. G. C. G. Eo 0 Eo 1, E. 1, E. 1; E. O, E. 1, 1 A. No OR:. Co G. 1,TRANSFER TO S o AND. B. E. 0, QUAD.(2.*MU,1,N), QUAD. (2.*MU, -4,N-l) QUAD,(2o*MU,-2,0) QUADo (LAMBDA,, 10) I QUAD. (LAMBDA,3, 5,N-) S FUNCTION RETURN FLOATING POINT EPSI, MU, LAMBDA NORMAL MODE IS INTEGER PROGRAM COMMON DUMMY (8004), PCH, PRSW, MXITER, EPSI, LD.ELTA, INIT, RUN, TSW, NAME(1),MAXJ VECTOR VALUES NAME = $EXAMPLES$ INTERNAL FUNCTION STATE.(X) = A + B + 5* C + 1 END OF FUNCTION 84

SA4SSFMBLE tEXEC-U-Tf-P PUNCH O3BJECTI,-D MP.1/0 DUMP. QUAD (OO0 _-_0 225_. _ 39/20/65 9 08 8.6 PM PAGE 1 (MOD 05 MAY 1965) UNIVERSITY OF MICHIGAN 20 SEP 1965,PRO6RAM LENGTH (OCTAL) 00352 wXrs.IWT Al lrfk.AA-CC],.CirTnO AAAn on PROGRAM COMMON LENGTH 17520 TrAn: A KI C r- r-O V 1 0rC - _? T. %. % uA N LOWEST ERASABLE DEFINED 77626 In QUAD 000 ti11T.tI'L IKAnIrItK ytuLIU UuuVV __IIKAN3t-K Vt'. IUK LFNG I.H UUUU_ Z _______ _.___ PROGRAN ENTRY POINT LOCATIONS AND NAMES 00002 SETUP 00031. U__ _ AD_ _ __ 00105 DIAG EXTERNAL ROUTINES CALLED BY THIS PROGRAM 00000 334751314563 00..PRIN'T 00001 255151465160 00 ERROR PROGRAN LISTING (NUMBERS IN OCTAL),'E__.___.] _...... ___.____. _.-..__..__.__ ___._ENRTQY__ SETIUPGQUA D D AG_ 00002 0634 00 1 00025 010 SETUP SXA XRA,1 00003 0534 00 1 77631 010 LXA Z,1 00004 0500 00 1 00250 010 CLA ADR,1 _00005 0621 00 00023 0.O._ __ L.. _ 00006 0621 00 0 00075 010 STA VB 00007 0621 0 _. 00076 010I S A ST__ 00010 0621 00 0 00122 010 STA COND 00011 0621 00 0 00131 010 STA VD 00012 0621 00 0 00132 010 STA VE 00013 621_(-a ___ 10 STA..__. _ HIGH _ ______ 00014 0621 00 0 00142 010 STA VF 00015 _4Q 00_0_035_. _ ADD_1 _ = 00016 0621 00 0 00160 010 STA VG 00017 0774 00 1 00001 00 AXT 1,1 00020 0634 00 1 13721 00 SXA TRANS,1 000221_ 0634 00 1 00251.910 _ X_ SX- _ SW, _ ____ 00022 0774 00 1 11610 00 AXT 5000,1 00O21.O_Ql_ _____ O ___OYA __ Z _ * _ __1 00024 2 00001 1 00023 010 TIX VA,1,1 00025 0774% 00 1 00000 00 XRA AXT **,1 00026 0500 60 4 00COi 00 CLA* 1,4 -..Q. 0027 __._0_1_ [.__Q0.4._0_Q. _ _____ STA FN ____ 00030 0020 00 4 00001 00 TRA 1,4 00031 0500 00 0 00251 010 QUAD__ C.._LA __ W _ __ 00032 0402 00 0 00351 010 SUB =1 00033 -0100 00 0 00171 010 TNZ ERSETA 00034 0634 00 4 00045 010 SXA XRB,4 0035__3_4__Q Q_ _ 2..0.POO_______ SXA _ XRC2_____ 00036 0634 00 1 00103 010 SXA XRD,1 0003_ __05_34 00._ 21371____.... _ _._ _LXA_ _TRANS 2____ __ 00040 3 03720 2 00201 010 TXH ERC,2,2000 00041 0500 60 4 33031 00 CLA* 1,4 00042 0601 00 2 13721 00 STO TRANS,2....0.0043 __ 63_4.O_2_ 034.6__010__.. SA _ _ _XRR+1_. ___ 00044 0074 00 4 00000 00 FN TSX **,4..0__45 0174 00 4 O0000_ X AXT. _ ** 4_ _ _ __ 00046 0774 00 2 00000 00 AXT **,2 00047 0340 00 0 00350 oln CAS =5000( 00050 0020 00 0 00212 010 TRA ERD 0.. 0051. _20.95...... _ _ _. __ _I_:RAK ___ 2......_ 00052 -0120 30 0 00220 0i10 TMI ERE ID QUAD 001 ID QUAD 002 ID QUAD 003

PAGBE 2 _UNIVERSITY OF MICHIGAI. 20 SEP 1965 00053 0601 00 2 17642 00 STO FIRST,2 00054 0400 60 4 00002 00 ADD* 2,4 00055 0601 00 2 23563 00 STO DEST.2 00056 0402 60 4 00002 00 SUB* 2,4 00057 0601 0 0 7 77717 _010 STO TMP 00060 0560 60 4 00003 00 LDQ* 3,4 00061 0200 00 0 27511 00 MPY DELTA 00062 0131 00 0 00000 00 XCA 00063 0400 00 0 77777 010_ ADD TMP 00064 0601 00 2 27504 00 STO LAST,2 00065 0734 00 _1 00000__00 00_ __ PAX 11 00066 -0634 00 1 00073 010 SXD THRU,1 00067 0534 O0 1 27511 00 LXA DELTA,1 00070 -0634 00 1 00077 010 SXD INC,1 00071 0500 00 2 17642 00 CLA _ F IRST 2 00072 0734 30 1 00000 00 PAX,1 00073 3 0003_0 1 00100 010 THRU TXH,___ Q 1.-** 00074 0502 00 2 13721 00 CLS TRANS,2 00075 0300 00 1 00000 30 VB FAD, *,1 00076 0601 00 1 00090 00 VC STO **,1 0P0017.. _1_0 1 _'O1 00073 010 INC TXI _ THRU,1,** 00100 1 00001 2 00101 310 BQ TXI *+1,2,1 00101 0634_00 2 13721_ 0 SXA TRANS,2 00102 0774 00 2 00000 00 XRC AXT **,2 00103 0774 00 1 00000 00 XRD AXT *i*,1 00104 0020 0.0 4 00001 00 TRA 1,4..Q5 5......0 00251__ D_ I.AG C.LA __SW _ 00106 0402 00 0 30351 310 SUB =1 0..l..7.-OlQQ-0 0. 00 _ 0 0_...0... ____ TZ _ ER SETB 00110 0634 00 1 -3164 010 SXA XRE,1 00111 0634 00 2 00165 310 SXA XRF,2 00112 0634 00 4 00166 010 SXA XRG,4 0._01.l3_..07. 0 0.1.000 O ____ AXT 1,1 00114 0534 00 4 27511 00 LXA DELTA,4 O.l _5. __-634 0 4 03.127. 10_.. SXD BKD,4L __ 00116 -0634 00 4 00156 010 SXD BKZ,4 00117 0534 00 4 13721 '0 LXA TRANS,4 ID QUAD 004 0\ 0'~ 00120 3 33720 4 00201 010 TXH ERC;4,2000 __001.21 _..22_..1 3 __62.LD......0............STK..... K+1___ 00122 -0520 00 1 3OOOG 00 COND NZT **,1 0012.3.1 0001 1.00122.10. TXI.*-1. 1. 00124 0754 00 I 09000 00 PXA,1 00125 0734 00 2 000OO 030 PAX,2 00126 3 11610 1 00164 010 TXH XRE,1,5000 00127 1 00000 2 00130 010 BKD TXI *+1_2. 00130 3 00000 2 00140 010 TXH HIGH,2,** 00131. 0500 00 1 C33C5 0 0'. 3 VD C.LA......__ 00132 0302 30 2 00000 00 VE FSB **P,2 001_3..Q7-6 QP.. _... ff._.Q 3:3__.Q. P_______ p ID QUAD 005 00134 00135 00136 00137 00140....014. _. 00142 00143 00144 0340 300 0020 30 3 0023 00 3 0020 00.0 0500 00 1.0) i. 0.. 4 0600 00 1 0754 00 1 0691 DO 4 27513 00140 ) 1 27.00127 00000 13721 310 010 30;^ CAS TRA TRA TRA HIGH CLA STO EPSI HIGH PKD BKD TRANS4 TRaLNS.4 02000 00 C0030 00 17642 00 VF STZ PXA STO **,1 *1 FI RST 4 ID QUAD 006

PAGE 3 00145 0601 00 4 23563' 00 00146 0754 00 2 00000 00 00147 0402 00 0 27511 0_ 00150 0601 00 4 27504 00 00151 1 00001 4 0f152 010 00152 3 03720 4 00201 010 00153 1 00001 1 00154 '010 00154 -0634 00 2 00157 010 OOl55 0634 00 1 00162 Oni 00156 1 00000 1 00157 010 00157 3 00000 1 00162 010 00160 0600 00 1 00000 00 00161 0020 00 0 O156 010 00162 0774 00 1 00000.00 00Q63 0020 0 0 00122 010 00164, 0774 00 1 00000 00 _00d 6 0774 00 2 000 0Q ---_ 00166 0774 00 4 00000 00 00167 0600 00 0 00251 010 00170 0020 00 4 00001 00 00171 0074 00 4 o00o o010 00172 -1 00000 0 00252 010 0073 -1 0000 0 _COOCOQ 00 00174 0074 00 4 00001 010 00175 on74 00 4 00000 010 00176 -1 00000 0 00261 010 00177 -1 00000_0 00000 00 00200 0074 00 4 00001 010 _00201 _74 _.. l. 3720 00 00202 0074 00 4 00000 010 00203 -1 00000 0 00270 010 00204 -1 00000 1 13721 00 _00205_-1 00000 1 17642 00 00206 -1 00000 1 23563 O0 _00207_-1 _I.000 27504 00 00210 -1 0000 0 00000 00 00211 0074 00 4 00001 010 00212 0601 00 0 77777 010 00213. 004.._00 4.._00 00 10 00214 -1 00000 0 00306 010 OQQ0215.._ — _..0QQ0. _ 00Q-_77Z71I- Q 00216 -1 00000 3 00003 00 00217 0020 00 0 00725 010 00220 0601 00 0 77777 010 _ D0221..074... 0 4...0 4 000....1.. 00222 -1 00000 0 00315 010 oQ2__. -L OQ0 O._. 0_777 010.... 00224 -1 00000 30 0000 00 _0225 0534 30 1 13721 30 00226 -3 00001 1 00246 010 00227 __1 777.. L 2..23. la_ 00230 -0634 00 1 00235 010 _.002 31 -...- O7.__ 0 _._0.....Q. 10. 00232 -1 00000 0 00324 010 __00233 0630 00 0 77777 010 00234 0774 00 1 00001 00 00235.3...OD0 1.00_ 245....0._ 00236 -1 00000 1 13721 00 UNIVERSITY_ QF MICHIGAN 20 SEP 1965 STO PXA.__. SUB,STO ___ TXI TXH TXI SXD SxXA.__ BKZ TXI STPZ.__TXH _ VG STZ TRA ENDZ AXT ___R__ XRE AXT X-RF. AX. XRG AXT STZ TRA FRSFTA PR I IT nDST.4,2 DELTA LAST,4 * l A+L. _ERC,4.2000 *+1,,]9 STPZ,2 ENDZ; L__. STPZ.t, 1,** EMnDZ 1.~** **,1 BKZ *, 1 *:* 1 -**.2 ** 4 SW "14 FA. CALL ERROR ERSETB PRINT FB,O CALL ERROR ERC AXI... _.. 2000,91 PRINT FC TOP TRANS,1 ___ _IlP ___ LRSISTL__ IOP DEST,1 ____ P_ LAST-_L__4 ENDIO CALL ERROR ERD STO TMP _PR I NT FD, ITM P.O_ TRA T8L ERE STO TMP -_ ___PRIN._L__FE, TMP.,___. TBL LXA TRANS,1 TXL END+1,l,1 _ _ I_ TXI ___* + 1_!_1_-1_ _ _.. SXD MAX,1 _...__..PRX__... P_ FTRL__ __ STZ TMP ID QUAD 007 ID QUAD 008 AXT 1,1 -_.IXH _..._..__END i, —** IOP TRANS,1

O0 0C PAGE 4 UNIVERSITY 00237 -1 00000 1 17642 00 IOP FIRST,1 00240 -1 00000 1 23563 00 IOP DEST, 1 -00241 -1 00000 1 27504 00 IOP LAST,1 -00242 0634 00 1 7?777 010 SXA TMP,1 00243 -1 00000 0 77777 010 lOP TMP 00?44 1 00001 1 00235 010 TXI MAX,l,l 00245 -1 00000 0 00000 00 END ENDIO 00246 0074 00 4 00001 010 CALL ERROR 00247 00000 0 66016 010' PZE V-5000 00250 0 0000 0 77626 010 ADR PZE V 00251 0 00000 0 00000 00 SW PZE 00252 030430105454 00 FA BCI 7,34H8**u* SETUP NOT CALLED BEFORE QUAD* 00253 545460622563 00 00254' 644760454663 00 00255 602321434325 00 ____ 00256 246022252646 00 00257 512560506421 00 0_ _ 00260 245460606060 00 00261 030430105454 00 FB BCI 7,34H8**** SETUP NOT CALLED BEFORE DIAG* 00262 545460622563 00 00263 644760454663 00 __ -.00264 602321434325 00 00265 246022252646 00 _ 00266 512560243121 00 00267 275460606060 00 00270 030130105454 00 FC BCI 9,31H8**** TRANSITION TABLE EXCEEDED/28HO**** LAST QUADR 00271 545460635121 00 _ __ 00272 456231633146 00 00273 456063212243 00 00274 256025672325 00 00275 252425246102 00 00276 103000545454 00 00277 546043216263 00_____ 00300 605064212451 00 00301 644743256066 00 __ BCI 5.UPLE WAS.../F45.5r3(1H.I19)* 00302 216260333333 00 00303 612604053305 O00 00304 730374013073 00 Q.305 _ 310111345460 O_. ____ __ 00306 011030105454 00 FO BCI 7,18H8**** STATE NUMBERI12,13H EXCEEDS 5000* 00307 545460626321 00__ _ _ _ _ _ ___ ___ 00310 632560456444 00 00311 222551310102 00 00312 730103306325 00. 0313 _..3Z 2232 5_2524_6Z. _____..... _._._ ___..__._ _.__ _ 00314 600500000054 00 40Q315.......Q.U30_0105aI 4 ___ _ _. FE _..L____ _C 7.,18H8**** STATE NUMBER_12,12H IS NEGATIVE _ 00316 545460626321 00 00317 632560456444 00 00320 222551310102 00 00321..... 73_0.1023Q031.... _ - b l 0__ _ _ _ _ _.__ 00322 626045252721 00 _-00323. 633 165Z55 6Q_..... __. _ _________ _ _ 00324 020530012364 0o FTBL RCI 9,25H1CURRENT TRANSITION TAPLE/1HOS39,5HTRANSSl5"5HFIRST 00Q25 5152545636 OQ 00326 635121456231 30 -_Q327 633_146456063 13' 00330 212243256101 30 ID QUA -009 ID QUAD 010 ID QUAD 0l OF MICHIGAN 20 SEP 1965

.PA1GE._S. __....-.._..._..... -..... _..... -..... -. -.. UNIVERSITY OF MICHIGAN! 20 SEP 1965 _00331 300 Q0_Z6I _ L7_3 00332 053063512145 0..0..... 33.. 26201 573.05.. 00334 302631516263 -..-_ 003.35 _._-6231.l7 3 043.... 00336 242562636201 00337 067304304321 00340 626362010373...... - 0341.0.10630.50 6421. 00342 245164474325....0343.. 604564442225 00344 516174013000 nnAr, 72AnsLr4r'All nq' __ O0 -_____ - X I o - ^- -- J - J. 00346 730431020034......... 03.4.7.......5_460..6 C6... 777-.............................. 776: 7765 _-___- __-__137Z 176z 23 5. 275(.......... 275: 275j,______... 5 PROGRAM LITERALS 00350. 000000011610 00351 0000300000001 *.** 00352 IS THE FIRST Lr 00..Q...g...... _..........................9............................................ 00 00........ Z..................... 00 00 BCI 2,,4I20)*...... 0.............................................. 77 TMP ERAS l101 31......... Z........ E R A S....2....... 26 V ER.AS 21_._-_.. _ _-RA.IR _P G4tml2......l Q.01 42 FIRST PGMCOM 2030 63. DEST PGMCOM.2000 04 LAST PGMCOM 2000 10. EP SI PGMCOM..3 11 DELTA PGMCOM LZ.__ _._M.AXJ__._PG.M _G____ __ END.NUMER/ (1HC F45. 5 ID QUAD 012 00 \O 'O O............ 0 3CATION NOT USED BY THIS PROGRAM. REFERENCES TO SYMBOLS DEFINED IN PROGRAM. TYPE SYMBOL VALUE REFERENCES TO SYMBOL (SYMBOL TYPES — R - RELOCATABLE, A - ABSOLUTE, F - ERASAPLE, T - T'iANSFER VECTOR) R ADR 00250 0,000.j4,253 -R._BKD a_.i_. _.C..L OIl 12.!~1.27..0.L36,,01n_37 _ _ R BKZ 00156 001160 1) 56 00161 R. BQ. 00100 00 73,0 1. R COND 00122 C0010, r 0122 00163 A DELTA 27511 00C61, 3' 067,,0114,00147,,00-35 A DEST 23563 00055, 3145,5 02C6,:C024) C,003503 R DIAG 00105 _.0C., r:2._p5....__ R END 00245 0226,:.)235, 0245 R ENDZ 03162 "0155 C'157 C. 162 A EPSI 275107 0134,') 50 R ERC 00201 C0r40, t) 12C,.00152,,CJ201 R ERD 00212 0.'35,1?0212 R ERE 00220'_ C 5r," 2 '2 T ERROR 300001 3' 174,,3 3321 1,0.,246 R ERSETA 00171 00-733 'I7171 R ERSETB 00175 'i"j7, '7175 R FA 00252 0j172,j252 R FB 03261 l0176, '0).261 R.....FC. _.._..2. 7 )O...(... 3.:?:.7 _.. R FD 00306 00 14, )3 06 R FE 00315 rC0222, 9315 A FIRST J7' +2 009;-53, '7 1, 001 44,.: 3235,.00237. 0'35 '

PAGE 6 UNIVERSITY OF MICHIGAN 20 SEP 196 R FN 00044 00027,'00044 R FT6L 00324 00232,00324 R H[GH 00140 00013,.00130,00135,00140 R -INC 00077 00070,00077 I.PRINT 00000 00171OO0175.00202.t00213,00221,00231 A LA.T 27504 00064,00150,00207t.00241 t00350 R I^X 00235 00230,.00235,00244 A IAX/J 27517 00350 R QUAD 00031 00000,.00002.00031 R SETUP 00002 00000,00002,00002 R STPZ 00157 00154,00156,00157 R SW 00251 00021,00031,00105,00167,.0025L R TOt 00225 00217,.00225 R THRU 00073 00066,00073,00077 E TUP 77777 00057,00063,00212,00215,00220,,00223,00233.00242,00243, 00350 A TRANS 13721 00020, 00037,00042, 00074,.00101,.00117,00141,.00204,00225,00236,00350 E V 77626 00247,00250 00350 R VA 00023 00005,00023.,00024 R Vs 00075 00006.,00075 R VC 00076 00007,00076 R VD 00131 00011'00131 R VE 00132 00012,00132 R Vs 00142 00014*,00142 R VS 00160 00016,.00160 R KRA 00025 00002:00025 R XRB 00045 00034,00043,00045 ft R XAR 00102 OQ35.00102 o R XR~ 00103 00036',00103 R XRE 00164,00110P,00126,00164 R XR6 00165 00111,00165 R XKR ool66 00112.00166 E Z 77631 00003,.00350 REFERENCES TO PROGRAM LITERALS eCE REFERENCES 00350 00047 00351 00015.00032i00106 FLAFBGE LiNES FATAL 0 NON-FATAL 0 **** DEFIN-IT-IONS' SYMBOL 55 OPERATION 0 MACRO 0 RMT 0 MACRO GALLS CORRECT 0 INCORRECT 0 **** CARD COUNTS INPUT 177 OUTPUT 13 5

UNIVERSITY OF MICHIGAN IIIIi 1'1111111111111Hlllllll 3 9015 03627 8284