dispatched to r servers, each having its own queue, according to the dispatching discipline (e.g. to the first r servers with least residual workload). The r servers act independently and serve the results according to some queueing discipline (e.g. First-Come-First-Served). After being processed, each result is routed to a single unlimited join queue. When the quorum q of results of the same workunit is obtained (and waits in the join queue), the workunit is marked as completed, and the unused results, which are either waiting in the queues, or being served by corresponding servers, are abandoned. Such models, known as (n, r, q) Fork-Join systems, have been studied extensively in [Jos16]. Note however, that, despite the simplicity of the model, only few analytic results are available and that too under stringent conditions (e.g. for the (n, n, n) M/M-type Fork-Join system [Jos16]), and in most cases only the asymptotic bounds or approximations are obtained [BDVD98,BM97]. In particular, the upper bound on the performance of the (n, r, r) Fork-Join system in [JSW15b] is done with the help of another closely related type of model, known as Split-Merge, or Split and Match queues. The key property of the Split-Merge system is that the service of the results of a particular workunit is started only after the previous workunit is completed (note, that in particular if q = 1 and r divides n, then Split-Merge and Fork-Join models coincide). The Split-Merge queues have been extensively studied since three decades ago. Two-server Split-Merge system was considered in [RP85], the departure process from Split-Merge was studied in [Rao90]. An emerging interest to the Split-Merge systems is related to new applications in Cloud and distributed systems [FL15].

In this work we suggest a Split-Merge model of the replication and quorum mechanism of a DG, which we refer below as (n, r, q) Split-Merge model. Our work extends the (n, n, n) Split-Merge model studied in [FL15], and elaborates the Fork-Join type models studied in [Jos16]. We give exact solutions for particular cases of the model that are analytically tractable, and study the effect of parameters, r and q, on the performance of the model by means of simulation.

The work is organized as follows. In section 2 we give a brief description of the DG technology, providing necessary details for the BOINC-based DG (the standard software used for distributed computing), and discuss the limitations of our approach. In section 3 we present the model, providing the necessary notation. In section 4 we validate the model and discuss the results of numerical experiments. The conclusion and possible extensions of the model are given in section 5.

DG technology is a distributed computing technology, that allows to utilize the idle resources (central processor unit time, memory and disk space) of a host in order to complete the computational project that consists of loosely coupled workunits. The hosts are either provided by volunteers (the Volunteer Computing, VC [NDS14]), or are a part of some controlled environment, e.g., the computational resources of a company (the Enterprise Desktop Grid, EDG [Iva15]). The project management, workunit generation, replication and result validation via quorum mechanism is performed by the management server of the project. During the project evaluation, the server produces workunits, then generates the necessary number of replicas of a workunit, and dispatches them to the hosts, either on their requests (VC), or impassively (EDG). The server receives and validates the results of computation, and assimilates them. However, it also has to deal with several key difficulties of VC: host availability, trust, malicious activity etc. [NDS14]. We also note, that the VC project requires a lot of social work with the volunteers, however, discussion of these issues is beyond the scope of this paper.

The EDG is a far less studied concept, which, however, has several simplifications compared to the VC. Namely, in most cases the hosts of an EDG may be treated as trusted, reliable, controllable and available. Moreover, since the hosts are directly under control of, say, the enterprise owner, the workunits may be PUSHed to the hosts (i.e. immediately dispatched), rather than PULLed by them (i.e. requested by clients in asynchronous regime, the basic technology used in the VC), and moreover may be canceled on demand [YJKA11]. In the present paper we mostly consider the model of EDG performance controlled by replication and quorum. However, some results of the study may be related to the VC model as well.

We briefly describe the key performance metrics, which are basically studied for the DG [ETA09]:

throughput -the number of workunits completed per unit time;

latency -the time from workunit generation (arrival) to workunit completion (departure);

delay -the time from workunit generation to the beginning of results computation;

starvation -the fraction of the idle time of the hosts, which is not used for computation; completion time -the average time it takes to complete a fixed number of workunits.

Moreover, the centralized structure of the DG network in some cases makes the project management server become the bottleneck of the system. Indeed, all the tasks related to dispatching, scheduling, result validation and assimilation, are done on a single server. In the current paper we mainly focus on the throughput, latency and delay of the DG model.

We consider an EDG project, that has a fixed number, say n, of identical hosts processing the workunits, generated by management server at instants of a point process with intensity λ. The workunit is immediately replicated and dispatched upon arrival to the hosts that have the least work left to be finished. Let r n be the number of results corresponding to a single workunit requested by the management server from the hosts, and let q r results be required to constitute a quorum. We assume, that the hosts start evaluating the results of a single workunit simultaneously, and the service time of results on each host is iid with distribution function F S . Immediately when q results are obtained, computation of r − q incomplete results (if any) is canceled by the management server. We assume that once received, the result is correct, i.e., the earliest q results received are valid. We also assume an unbounded time for result computation, that is, q results are eventually received. We are interested in the performance metrics of such a computation.

Let S 1 , . . . , S r be iid random variables corresponding to the (potential) times of computation of the results of a single workunit. Since the results are started simultaneously, then the time to complete evaluation of the workunit is distributed as the q th order statistics of r r.v., which we denote by S q:r . Recall, that

where F S (x) := 1 − F S (x). In particular, the minimum of S 1 , . . . , S r is distributed as

Thus, the system under study is equivalent to an M/G/ n r queueing system with general service time S q:r . Note however, that it is not necessary, that F Sq:r (x) := P (S q:r x) belongs to the same class of distributions, as F S . Obtaining the distribution of q th order statistics analytically is difficult in general. Note that if the service times are exponentially distributed with F S (x) = 1 − e −µx , µ > 0, then S q:r is exponentially distributed only for q = 1, which provides F Sq:r (x) = 1 − e −rµx , r 1.

The class of distributions, that is closed under the operation of obtaining order statistics, is phase-type (PH) distribution. A continuous-time PH distribution of order m with representation (β, B) is the time to absorption of the continuous-time Markov chain with a fixed number m of transient states (phases), initial distribution of states β = (β 1 , . . . , β m ) and transition subintensity matrix B. The vector b 0 = −B1 0 gives the intensity of transitions to absorbing state m + 1 (for more details on this type of distributions see e.g. [He14,Neu81]). Recall also, that if S has a (β, B) PH distribution, then ES = −βB −1 1 and ES 2 = 2βB −2 1 [BKF14]. Let F S be a (β, B) PH distribution. Now we give a procedure to obtain PH representation of the distribution F Sq:r (x) (a detailed discussion on the procedure in case of non-identical PH distributions may be found in [BN17]).

First, we need to extend the phase space. Denote M k := {1, . . . , m} k the set of m k lexicographically ordered k-tuples corresponding to the states of k processes that have not yet reached absorbing state. Now we construct the process

that corresponds to evolution of the r processes until q absorptions, where the number of components equals r + (r − 1) + • • • + (r − q + 1). Below we explain the evolution of the process in detail. Initially the first r components of the process are independent copies of r.v. with PH(β, B) distribution. Thus, the evolution of the process (X 1 (t), . . . , X r (t)) ∈ M r without absorptions is governed by the m r × m r subintensity matrix B ⊕r , where

Intuitively the notation above means, that only one of the r.v. X 1 , . . . , X r makes its own independent transition, that does not end into absorption. When an absorption occurs, we exclude the absorbing component, switching to the process (X r+1 , . . . , X 2r−1 ) ∈ M r−1 . Thus, the transitions of the process related to absorption of one of the components are governed by m r × m r−1 matrix B (r) 0 . We stress, that the rate [B (r) 0 ] i,i of transition from a state i = (i 1 , . . . , i r ) ∈ M r to the state i = (i 1 , . . . , i r−1 ) ∈ M r−1 depends on the number of possibilities to remove exactly one component of i and arrive at i as follows

where N (i, i ) = {t r − k + 1 : i j = i j , 1 j t, i j+1 = i j , j > t} is the set of indices of such components of i, that, being removed from i, convert it to i . We note, however, that, by exploiting the lexicographical order of M r and M r−1 , a more straightforward expression for the transition rate matrix appears

(with obvious conventions if r < 2) where I j is the identity matrix of size j 1.

Then the transitions of the process X (t) are governed by the following bidiagonal matrix (where zeroes are the zero matrices of the corresponding dimension)

It remains to note, that the initial distribution of the process is given by

and the transition intensity to the super-absorbing state corresponding to exactly q absorptions of r independent processes is given by −B q:r 1. It may be seen, that row sums of the matrix B q:r are zero for all rows, except the last m r−q+1 , which correspond to the phases of process with q − 1 absorptions. It may be now seen, that the r.v. S q:r has a PH(β q:r , B q:r ) distribution.

The PH representation allows to obtain immediately the stability condition and performance measures of the system as follows.

Stability condition of the system

where, recall, ES q:r = −β q:r B −1 q:r 1. Thus, the maximal throughput of the system equals

The explicit forms for the mean stationary delay (ED), as well as the mean latency (EW ), are available for r = n (for more details see Appendix E in [Jos16], case q = r = n is considered in [FL15]), since the system in this case is equivalent to an M/P H/1 system:

For r < n it is possible either to use approximation [Jos16], or rely on simulation. Alternatively, the solution of continuous-time QBD process corresponding to the M/P H/ n r system may be obtained in terms of the intensity matrix R, which (except some special cases) is derived only numerically. However, discussion of this approach is beyond the scope of this paper and will be presented elsewhere.

We note, that it is easy to sample exactly from q th order statistics of r independent PH(β, B)-type variables by sampling from a single PH(β q:r , B q:r ) r.v. The drawback of this approach is the growing size of the number of phases that requires to process square matrices of order q−1 i=0 m r−i = (m r+1 − m r−q+1 )/(m − 1). Note however, that it is mostly a technical limitation. Now we focus on the latency in some particular cases, that are most analytically tractable. An important particular case is the M/M/ n r -type system with exponential interarrival times. This corresponds to an (n, r, 1)-type replication with exponentially distributed service times (with intensity µ). In this case we have (n, r, 1) split-merge model, i.e., we send r results and wait for the fastest, we obtain the Exp(rµ) service distribution of the minimum of r exponentials, and obtain the M/M/ n r -type model. That is, the sojourn time equals

and in particular, if r = 1 (no replication at all), the sojourn time equals

where

is the Erlang C-formula. On the other hand, if r = n (full replication), we obtain the M/M/1 system with service intensity nµ, and the mean stationary sojourn time is EW = 1 nµ + 1 nµ−λ .

In this paper we have presented a model of the EDG system with (n, r, q)-type replication and quorum mechanism. The results of numerical experiments illustrate the practical applicability of the model and show good correspondence to the obtained, as well as earlier known, theoretical results. Note however, that the numerical experiments mostly cover the so-called log-convex distributions. It is known, that the monotonicity of dependence of the mean stationary latency on the values q, r may the opposite, compared to the log-concave distribution (for a detailed discussion on the impact of log-concavity/convexity on latency for the Fork-Join type models see [Jos16]). Moreover, the dependence may not be monotonic for the distributions with the so-called heavy tails (e.g. Pareto distribution). Thus, a more intensive numerical study is required, and we leave this for future research.

We also note, that an extension to heterogeneous servers may be done taking into accounts results of [AM01], and order statistics of heterogeneous PH distributions are discussed in [BN17]. 6) for parameter sweep experiment, with r = 1, . . . , 5 and q = 1, . . . , r.

[BM97] Simonetta Balsamo and Ivan Mura. On queue length moments in fork and join queuing networks with general service times, pages 218-231. Springer Berlin Heidelberg, Berlin, Heidelberg, 1997. q q q q q q q q q q 2 4 6 8 10 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Replication Latency q q q q q q q q q q q q q q q q q q=1 q=2 q=3 theoretical

Figure 3: Comparison of the sample mean latency for q = 1, 2, 3, r = 1, . . . , 10 and n = 100 as a result of parameter sweep experiment with M/G/ n r queueing system with general service times S q:r observed serving 5 • 10 5 tasks, compared to the theoretical value of average stationary latency EW obtained from (7).