A linear parallel algorithm to compute bisimulation and relational coarsest partitions

The most efficient way to calculate strong bisimilarity is by calculation the relational coarsest partition on a transition system. We provide the first linear time algorithm to calculate strong bisimulation using parallel random access machines (PRAMs). More precisely, with $n$ states, $m$ transitions and $|\mathit{Act}|\leq m$ action labels, we provide an algorithm on $max(n,m)$ processors that calculates strong bisimulation in time $O(n+|\mathit{Act}|)$ and space $O(n+m)$. The best-known PRAM algorithm has time complexity $O(n\log n)$ on a smaller number of processors making it less suitable for massive parallel devices such as GPUs. An implementation on a GPU shows that the linear time-bound is achievable on contemporary hardware.


Introduction
The notion of bisimilarity for Kripke structures and Labeled Transition Systems (LTSs) is commonly used to define behavioural equivalence. Deciding this behavioural equivalence is important in the field of modelling and verifying concurrent systems [4,15]. Kanellakis and Smolka proposed a partition refinement algorithm for obtaining the bisimilarity relation for Kripke structures [11]. The proposed algorithm has a run time complexity of O(nm) where n is the number of states and m is the number of transitions of the input graph. Later, a more sophisticated refinement algorithm running in O(m log n) steps was proposed by Paige and Tarjan [16].
In recent years the increase in the speed of sequential chip design has stagnated due to a multitude of factors such as energy consumption and heat generation. In contrast, parallel devices such as graphics processing units (GPUs) keep increasing rapidly in computational power. In order to profit from the acceleration of these devices, we require algorithms with massive parallelism. The article "There's plenty of room at the Top: What will drive computer performance after Moore's law" by Leierson et al. [14] indicates that the advance in computational performance will come from software and algorithms that can employ hardware structures with a massive number of simple, parallel processors, such as GPUs. In this paper, we propose such an algorithm to decide bisimilarity.
Deciding bisimilarity is P -complete [1], which suggests that bisimilarity is an inherently sequential problem. This fact has not withheld the community from searching efficient parallel algorithms for deciding bisimilarity of Kripke structures. In particular, Lee and Rajasekaran [13] proposed a parallel algorithm based on the Paige Tarjan algorithm that works in O(n log n) time complexity using m log n log log n Concurrently Read and Concurrently Write (CRCW) processors.
In this work, we improve on the best known theoretical bound for PRAM algorithms using a higher degree of parallelism. The proposed algorithm improves the run time complexity to O(n) on max(m, n) processors and is based on the sequential algorithm of Kanellakis and Smolka [11]. The larger number of processors used in this algorithm favours the increasingly parallel design of contemporary and future hardware. In addition, the algorithm is optimal w.r.t. the sequential Kanellakis-Smolka algorithm, meaning that overall, it does not perform more work than its sequential counterpart.
We first present our algorithm on Kripke structures where transitions are unlabelled. However, as labelled transition systems (LTSs) are commonly used, and labels are not straightforward to incorporate in an efficient way (cf. for instance [21]), we discuss how our algorithm can be extended to take action labels into account. This leads to an algorithm with a run time complexity of O(n + |Act|), with Act the set of action labels.
Our algorithm has been designed for and can be analyzed with the CRCW PRAM model, following notations from [20]. This model is an extension of the normal RAM model, allowing multiple processors to work with shared memory. In the CRCW PRAM model, parallel algorithms can be described in a straightforward and elegant way. In reality, no device exists that completely adheres to this PRAM model, but with recent advancements, hardware gets better and better at approximating the model since the number of parallel threads keeps growing. We demonstrate this by translating the PRAM algorithm to GPU code. We straightforwardly implemented our algorithm in CUDA and experimented with an NVIDIA Titan RTX, showing that our algorithm performs mostly in line with what our PRAM algorithm predicts.
The paper is structured as follows: In Section 2, we recall the necessary preliminaries on the CRCW PRAM model and state the partition refinement problems this paper focuses on. In Section 3, we propose a parallel algorithm to compute bisimulation for Kripke structures, which is also called the Relational Coarsest Partition Problem (RCPP). In this section, we also prove the correctness of the algorithm and provide a complexity analysis. In Section 4, we discuss the details for an implementation with multiple action labels, thereby supporting LTSs, which forms the Bisimulation Coarsest Refinement Problem (BCRP). In Section 5 we discuss the results of the implementation and in Section 6 we address the usage of weaker PRAM models. Finally, in Section 7, we discuss related work.

The PRAM Model
The Parallel Random Access Machine (PRAM) is a natural extension of the normal Random Access Machine (RAM), where an arbitrary number of parallel programs can access the memory. Following the definitions of [20] we use a version of PRAM that is able to Concurrently Read and Concurrently Write (CRCW PRAM). It differs from the model introduced in [8] in which the PRAM model was only allowed to concurrently read from the same memory address, but concurrent writes (to the same address) could not happen. We call the model from [8] an Concurrent Read, Exclusive Write (CREW) PRAM model.
A CRCW PRAM consists of a sequence of numbered processors P 0 , P 1 , . . . . These processors have all the natural instructions of a normal RAM such as addition, subtraction, and conditional branching based on the equality and less-than operators. There is an infinite amount of common memory the processors have access to. The processors have instructions to read from and write to the common memory. In addition, a processor P i has an instruction to obtain its unique index i. A PRAM also has a function P : N → N which defines a bound on the number of processors given the size of the input.
All the processors have the same program and run synchronized in a single instruction, multiple data (SIMD) fashion. In other words, all processors execute the program in lockstep. Parallelism is achieved by distributing the data elements over the processors and having the processors apply the program instructions on 'their' data elements.
Initially, given input consisting of n data elements, the CRCW PRAM assumes that the input is stored in the first n registers of the common memory, and starts the first P (n) processors P 0 , P 1 , . . . , P P (n)−1 .
We need to define what the behaviour of the machine will be whenever a concurrent write happens. The way to handle this memory contention in concurrent writes is usually by assuming one of the following: • (Common) All processors try to write the same value and succeed, otherwise, the writes are not legal and fail; • (Arbitrary) Only one arbitrary attempt to write succeeds; • (Priority) Only the processor with the lowest index succeeds in writing.
The algorithm proposed in this paper works if we make either the arbitrary or the priority assumption. In Section 6 we explain how we can adapt it to work under the common assumption.
A parallel program for a PRAM is called optimal w.r.t. a sequential algorithm if the total work done by the program does not exceed the work done by the sequential algorithm. More precisely, if T is the parallel run time and P the number of processors used, then the algorithm is optimal w.r.t. a sequential algorithm running in S steps if P · T ∈ O(S).
The computational complexity of these models is well studied and there is a close relation between circuit complexity and the complexity of PRAM algorithms [20].

Strong Bisimulation
To formalise concurrent system behaviour, we use LTSs.
Definition 1 (Labeled Transition System). A Labeled Transition System (LTS) is a threetuple A = (S, Act, →) where S is a finite set of states, Act a finite set of action labels, and →⊆ S × Act × S the transition relation.
Let A = (S, Act, →) be an LTS. Then, for any two states s, t ∈ S and a ∈ Act, we write s a − → t iff (s, a, t) ∈ →.
Kripke structures differ from LTSs in the fact that the states are labelled as opposed to the transitions. In the current paper, for convenience, instead of using Kripke structures where appropriate, we reason about LTSs with a single action label, i.e., |Act| = 1. Computing the coarsest partition of such an LTS can be done in the same way as for Kripke structures, apart from the fact that in the latter case, a different initial partition is computed that is based on the state labels (see, for instance, [9]). Definition 2 (Strong bisimulation). On an LTS A = (S, Act, →) a relation R ⊆ S × S is called a strong bisimulation relation if and only if it is symmetric and for all s, t ∈ S with sRt and for all a ∈ Act with s a − → s , we have: Whenever we refer to bisimulation we mean strong bisimulation. Two states s, t ∈ S in an LTS A are called bisimilar, denoted by s t, iff there is some bisimulation relation R for A that relates s and t.
A partition π of a finite set of states S is a set of subsets that are pairwise disjoint and whose union is equal to S, i.e., B∈π B = S. Every element B ∈ π of this partition π is called a block.
We call partition π a refinement of π iff for every block B ∈ π there is a block B ∈ π such that B ⊆ B. We say a partition π of a finite set S induces the relation R = {(s, t) | ∃B ∈ π.s ∈ B ∧ t ∈ B}. This is an equivalence relation of which the blocks of π are the equivalence classes.
Given an LTS A = (S, Act, →) and two states s, t ∈ S we say that s reaches t with action a ∈ Act iff s a − → t. A state s reaches a set U ⊆ S with an action a iff there is a state t ∈ U such that s reaches t with action a. A set of states V ⊆ S is called stable under a set of states U ⊆ S iff for all actions a either all states in V reach U with a, or no state in V reaches U with a. A partition π is stable under a set of states U iff each block B ∈ π is stable under U . The partition π is called stable iff it is stable under all its own blocks B ∈ π. Fact 1. [16] Stability is inherited under refinement, i.e. given a partition π of S and a refinement π of π, then if π is stable under U ⊆ S, then π is also stable under U .
The main problem we focus on in this work is called the bisimulation refinement problem (BCRP). It is defined as follows: Input: An LTS M = (S, Act, →). Output: The partition π of S which is the coarsest partition, i.e., has the smallest number of blocks, that forms a bisimulation relation.
In a Kripke structure, the transition relation forms a single binary relation, since the transitions are unlabelled. This is also the case when an LTS has a single action label. In that case, the problem is called the Relational Coarsest Partition Problem (RCPP) [11,13,16]. This problem is defined as follows: Input: A set S, a binary relation →: S × S and an initial partition π 0 Output: The partition π which is the coarsest refinement of π 0 and which is a bisimulation relation.
It is known that BCRP is not significantly harder than RCPP as there are intuitive translations from LTSs to Kripke structures [7,Dfn. 4.1]. However, some non-trivial modifications can speed-up the algorithm for some cases, hence we discuss both problems separately. In Section 3, we discuss the basic parallel algorithm for RCPP, and in Section 4, we discuss the modifications required to efficiently solve the BCRP problem for LTSs with multiple action labels.

A Sequential Algorithm
In this section, we discuss a sequential algorithm based on one of Kanellakis and Smolka [11] for RCPP. This is the basic algorithm which we adapt to the parallel PRAM algorithm. The algorithm starts with an input partition π 0 and refines all blocks until a stable partition is reached. This stable partition will be the coarsest refinement that defines a bisimulation relation. The sequential algorithm, Algorithm 1, works as follows. Given are a set S, a relation →⊆ S × S, and an initial partition π 0 of S. Initially, we mark the partition as not necessarily stable under all blocks by putting these blocks in a set Unstable. In any iteration of the algorithm, if a block B of the current partition is not in Unstable, then the current partition is stable under B. If Unstable is empty, the partition is stable under all its blocks, and the partition represents the required bisimulation.
As long as some blocks are in Unstable (line 3), a single block B ∈ π is taken from this set (line 4) and we split the current partition such that it becomes stable under B. Therefore, we refer to this block as the splitter. The set S = {s ∈ S | ∃t ∈ B.s → t} is the reverse image of B (line 6). This set consists of all states that can reach B, and we use S to define our new blocks. All blocks B that have a non-empty intersection with S , i.e., B ∩ S = ∅, and are not a subset of S , i.e., B ∩ S = B (line 7), are split in the subset of states that reach S and the subset of states that do not reach S (lines 8-9). These two new blocks are added to the set of Unstable blocks (line 10). The number of states is finite, and blocks can be split only a finite number of times. Hence, blocks are only finitely often put in Unstable, and so the algorithm terminates.

The PRAM Algorithm
Next, we describe a PRAM algorithm to solve RCPP that is based on the sequential algorithm given in Algorithm 1.

Block representation
Given an LTS A = (S, Act, →) with |A| = 1 and |S| = n states, we assume that the states are labeled with unique indices 0, . . . , n − 1. A partition π in the PRAM algorithm is represented by assigning a block label from a set of block labels L B to every state. The number of blocks can never be larger than the number of states, hence, we use the indices of the states as block labels: L B = S. We exploit this in the PRAM algorithm to efficiently select a new block label whenever a new block is created. We select the block label of a new block by electing one of its states to be the leader of that block and using the index of that state as the block label. By doing so, we maintain an invariant that the leader of a block is also a member of the block.
In a partition π, whenever a block B ∈ π is split into two blocks B and B , the leader s of B which is part of B becomes the leader of B , and for B , a new state t ∈ B is elected to be the leader of this new block. Since the new leader is not part of any other block, the label of t is fresh with respect to the block labels that are used for the other blocks. This method of using state leaders to represent subsets was first proposed in [24,23].

Data structures
The common memory contains the following information: 6. The following is stored in lists of size n, for each potential block with block label i: As input, we assume that each state with index i has an input variable I i ∈ L B that is the initial block label. In other words, the values of the I i variables together encode π 0 . Using this input, the initial values of the block label block i variables are calculated to conform to our block representation with leaders. Furthermore in the initialization, unstable i = false for all i that are not used as block label, and true otherwise.
Step 1: Select current block:= B s 4 Step 2: Mark nodes s 1 , s 2 Step 3: Split B s 1 into B s 1 , B s 3 Figure 1: One iteration of Algorithm 2

The algorithm
We provide our first PRAM algorithm in Algorithm 2. The PRAM is started with max(m, n) processors. These processors are dually used for transitions and states. The algorithm performs initialisation (lines 1-6), after which each block has selected a new leader (lines 3-4), ensuring that the leader is one of its own states, and the initial blocks are set to unstable. Subsequently, the algorithm enters a single loop that can be explained in three separate parts.
Splitter selection (lines [8][9][10][11][12][13][14], executed by n processors. Every variable mark i is set to false. After this, every processor with index i will check unstable i . If block i is marked unstable the processor tries to write i in the variable C. If multiple write accesses to C happen concurrently in this iteration, then according to both the arbitrary and the priority PRAM model (see Section 2), only one process j will succeed in writing, setting C := j as splitter in this iteration.
Mark states (lines [15][16][17], executed by m processors. Every processor i is responsible for the transition s i − → t i and checks if t i (target i ) is in the current block C (line 15). If this is the case the processor writes true to mark source i where source i is s i . This mark now indicates that s i reaches block C.

Correctness
The block i list in the common memory at the start of iteration k defines a partition π k where states s ∈ S with equal block labels block i form the blocks: A run of the program produces a sequence π 0 , π 1 , . . . of partitions. Observe that partition π k is a refinement of every partition π 0 , π 1 , . . . , π k−1 , since blocks are only split and never merged.
A partition π induces a relation of which the blocks are the equivalence classes. For an input partition π 0 we call the relation induced by the coarsest refinement of π 0 that is a bisimulation relation π 0 .
We now prove that Algorithm 2 indeed solves RCPP. We first introduce Lemma 1 which is invariant throughout execution and expresses that states which are related by π 0 are never split into different blocks. This lemma implies that if a refinement forms a bisimulation relation, it is the coarsest. Lemma 1. Let S be the input set of states, →: S × S the input relation and π 0 the input partition. Let π 1 , π 2 , . . . be the sequence of partitions produced by Algorithm 2, then for all initial blocks B ∈ π 0 , states s, t ∈ B and iteration k ∈ N: Proof. This is proven by induction on k. In the base case, π 0 , this is true by default. Now assume for a particular k ∈ N that the property holds. We know that the partition π k+1 is obtained by splitting with respect to a block C ∈ π k . For two states s, t ∈ S with s π 0 t we know that s and t are in the same block in π k . In the case that both s and t do not reach C, then mark s = mark t = false. Since they were in the same block, they will be in the same block in π k+1 . Now consider the case that at least one of the states is able to reach C. Without loss of generality say that s is able to reach C. Then there is a transition s → s with s ∈ C. By Definition 2, there exists a t ∈ S such that t → t and s π 0 t . By the induction hypothesis we know that since s π 0 t , s and t must be in the same block in π k , i.e., t is in C. This witnesses that t is also able to reach C and we must have mark s = mark t = true. Since the states s and t are both marked and are in the same block in π k , they will remain in the same block in π k+1 . Lemma 2. Let S be the input set of states with →: S × S, L B = S the block labels, and π n the partition stored in the memory after termination of Algorithm 2. Then the relation induced by π n is a bisimulation relation.
Proof. Since the program finished, we know that for all block indices i ∈ L B we have unstable i = false. For a block index i ∈ L B , unstable i is set to false if the partition π k , after iteration k, is stable under the block with index i and set to true if it is split. So, by Fact 1, we know π n is stable under every block B, hence stable. Next, we prove that a stable partition is a bisimulation relation.
We show that the relation R induced by π n is a bisimulation relation. Assume states s, t ∈ S with sRt are in block B ∈ π n . Consider a transition s → s with s ∈ S. State s is in some block B ∈ π n , and since the partition is stable under block B , and s is able to reach B , by the definition of stability, we know that t is also able to reach B . Therefore, there must be a state t ∈ B such that t → t and s Rt . Finally, by the fact that R is an equivalence relation we know that R is also symmetric, therefore it is a bisimulation relation. Proof. By Lemma 2, the resulting partition is a bisimulation relation. Lemma 1 implies that it is the coarsest refinement which is a bisimulation.

Complexity analysis
Every step in the body of the while-loop can be executed in constant time. So the asymptotic complexity of this algorithm is given by the number of iterations. Proof. In iteration k ∈ N of the algorithm, let us call the total number of blocks N k ∈ N and the number of unstable blocks U k ∈ N. Initially, N 0 = U 0 = |π 0 |. In every iteration k, a number of blocks l k ∈ N is split, resulting in l k new blocks, so the new total number of blocks at the end of iteration k is N k+1 = N k + l k .
First the current block C in iteration k which was unstable is set to stable which causes the number of unstable blocks to decrease by one. In this iteration k the l k blocks B 1 , . . . , B l k are split, resulting in l k newly created blocks. These l k blocks are all unstable. A number of blocks l k ≤ l k of the blocks B 1 , . . . B l k , were stable and are set to unstable again. The block C which was set to stable is possibly one of these l k blocks which were stable and set to unstable again. The total number of unstable blocks at the end of iteration k is U k+1 = U k + l k + l k − 1.
For all k ∈ N, in iteration k we calculate the total number of blocks

and the total number of newly created blocks is
The number of unstable blocks is always positive, i.e., U k ≥ 0, and the total number of blocks can never be larger than the number of states, i.e., N k ≤ n, so the total number of iterations z is bounded by z ≤ 2N z − |π 0 | ≤ 2n − |π 0 |.

Bisimulation Coarsest Refinement Problem
In this section we extend our algorithm to the Bisimulation Coarsest Refinement Problem (BCRP), i.e., to LTSs with multiple action labels.
Solving BCRP can in principle be done by translating an LTS to a Kripke structure, for instance by using the method described in [18]. This translation introduces a new state for every transition, resulting in a Kripke structure with n + m states. If the number of transitions is significantly larger than the number of states, then the number of iterations of our algorithm increases undesirably.

The PRAM Algorithm
Instead of introducing more states, we introduce multiple marks per state, but in total we have no more than m marks. For each state s, we use a mark variable for each different outgoing action label relevant for s, i.e., for each a for which there is a transition s a − → s to some state s . Each state may have a different set of outgoing action labels and thus a different set of marks. Therefore, we first perform a preprocessing procedure in which we group together states that have the same set of outgoing action labels. This is valid, since two bisimilar states must have the same outgoing actions. That two states of the same block have the same set of action labels is then an invariant of the algorithm, since in the sequence of produced partitions, each partition is a refinement of the previous one. For the extended algorithm, we need to maintain extra information in addition to the information needed for Algorithm 2. For an input LTS A = (S, Act, − →) with n states and m transitions this is the following extra information: 1. Each action label has an index a ∈ {0, . . . , |Act| − 1}. (c) nr marks s , the number of marks this state has, thus the length of list mark s .

mark length:
The total length of all the mark s lists together, i.e., the sum of all the nr marks s . This allows us to reset all marks in constant time using mark length processors. This number is not larger than the number of transitions (mark length ≤ m).

5.
In a list of size n, we store for each state s ∈ S a variable split s . This indicates if the state will be split off from its block.
With this extra information, we can alter Algorithm 2 to work with labels. The new version is given in Algorithm 3. The changes involve the following: 1. Lines 7-9: Reset the mark list.
3. Line 17: When marking the transitions, we do this for the correct action label, using order i . Note the indexing into mark . It involves the offset for the state source i , and order i .

Lines 19-21:
We tag a state to be splitted off when it differs for any action label from the block leader.
5. Line 24: If a state was tagged to be splitted off in the previous step, it should split off from its leader.
6. Line 29: If any block was split, the partition may not be stable w.r.t. the splitter.

Preprocessing.
To use the above algorithm, we need to do two preprocessing steps. First, we need to partition the states w.r.t. their set of outgoing action labels. This can be done with an altered version of Algorithm 2. Instead of splitting on a block at line 15, we split on an action a ∈ A. We visit all transitions, and we mark the source if it has the same action label a. This can be found in Algorithm 4. After executing Algorithm 4, each block can split in two blocks: a block that contains states that have a as an outgoing action label and a block with states that do not have this outgoing action label. After doing this for all different action labels we end up with a partition of blocks, in which all states of a block have the same set of outgoing action labels, and each pair of states from different blocks have different sets of outgoing action labels. Using m processors, this partition can be constructed in O(|Act|) time.
For the second preprocessing step, we need to gather the extra information that is needed in Algorithm 3. Only a i is part of the input, the others need to be calculated. We start our preprocessing by sorting the transitions by (source i , a i ), which can be done in O(log m) time with m processors, for instance using a parallel merge sort [5]. In order to calculate order i and nr marks s , we first calculate action switch i for each transition i, which is done in Algorithm 5. See Figure 2 for an example. Now, order i can be calculated with a parallel segmented prefix sum [19] (also called a segmented scan) of action switch. A parallel segmented sum can be performed on action switch to calculate nr marks s , where we make sure to set nr marks s to 0, if state s has no outgoing transitions. Finally, off s , for the mark offsets, can be constructed as a list and calculated by applying a parallel prefix sum on nr marks s . The code in Algorithm 5 takes O(1) time on m processors, and a parallel segmented (prefix) sum takes O(log m) time [19].
In total the preprocessing takes O(|Act| + log m) time.

Complexity & Correctness
For Algorithm 3, we need to prove why it takes a linear number of steps to construct the final partition. This is subtle, as an iteration of the algorithm does not necessarily produce a stable block. in l k new blocks, meaning that N k+1 = N k + l k . All new l k blocks are unstable and a number l k ≤ l k of the old blocks that are split, were stable at the start of iteration k and now unstable. If l k = l k = 0 there are no blocks split and the current block C becomes stable. We indicate this with a variable c k : c k = 1 if l k = 0, and c k = 0, otherwise. The total number of iterations up to iteration k in which no block is split is given by k−1 i=0 c i . The number of iterations in which at least one block is split is given by k − k−1 i=0 c i . If in an iteration k at least one block is split, the total number of blocks at the end of iteration k is strictly higher than at the beginning, hence for all k ∈ N, i=0 c i is an upper bound for k. We derive an upper bound for the number of iterations in which no blocks are split using the total number of unstable blocks. In iteration k there are U k = k−1 i=0 (l i +l i )− k−1 i=0 c i +|π 0 | unstable blocks. Since the sum of newly created blocks k−1 i=0 (l i ) = N k − |π 0 | and l i ≤ l i , the number of unstable blocks U k is bounded by This gives the bound on the total number of iterations z ≤ 3N z − |π 0 | ≤ 3n − |π 0 |.
With the time for preprocessing this makes the total run time complexity O(n + |Act| + log m). Since the total number of transitions m is bounded by |Act| × n 2 , this simplifies to O(n + |Act|).
Concerning correctness, we need to address two things. Firstly, as argued above, we start with a different partition compared to Algorithm 2, but it is a valid choice since states with different outgoing labels can never be bisimilar. Secondly, although the partition may not become stable w.r.t. the splitter, this will eventually occur, and the algorithm will only stop once the partition is stable w.r.t. all blocks. Therefore, the algorithm will produce the coarsest bisimulation relation.

Experimental Results
In order to validate the proposed algorithm, we implemented Algorithm 3 from Section 4. The implementation targets graphics processing units (GPUs) since a GPU closely resembles a PRAM and supports a large amount of parallelism. The algorithm is implemented in CUDA version 11.1 with use of the Thrust library. 1 As input, we chose all benchmarks of the VLTS benchmark suite 2 for which the implementation produced a result within 10 minutes. The Benchmark name n m |Act| |Blocks| #It T pre T alg T total #It/n #It/ |Blocks| T total /n T alg /#It VLTS benchmarks are LTSs that have been derived from real concurrent system models. The experiments were run on an NVIDIA Titan RTX with 24 GB memory and 72 Streaming Multiprocessors, each supporting up to 1,024 threads in flight. Although this GPU supports 73,728 threads in flight, it is very common to launch a GPU program with one or even several orders of magnitude more threads, in particular to achieve load balancing between the Streaming Multiprocessors and to hide memory latencies. In fact, the performance of a GPU program usually relies on that many threads being launched.
Our implementation is purely a proof of concept, to show that our algorithm can be mapped to actual hardware and to understand how the algorithm scales with the number of states and transitions.
In the implementation, we have to make a few adjustments, since a GPU differs in some aspects from a PRAM. To make memory updates globally visible, we need to synchronize at certain points of Algorithm 3, otherwise the changes in the memory are not consistent. We do this by splitting up the algorithm in different kernels (functions that execute in parallel on a GPU) since after a kernel run all processors (threads) are synchronized.
To be precise, in Algorithm 3 we need to synchronize after: • Line 15: To make sure the mark and split lists are reset and the splitter (C) is the same for all threads.
• Line 18: To make sure every thread has the same view of the mark list.
• Line 21: To synchronize the mark list.
• Line 25: To make sure the next leader for states that split off (new leader block i ) is chosen consistently among threads.
We have chosen to allow race conditions in our implementation, for instance at Line 6 where multiple blocks can mark themselves as current (C). Strictly speaking, this is not safe in the CUDA programming model, but it does work for 32 bit words. This can be easily adjusted using atomic instructions, although this will result in sequentializing write accesses to the same memory location, meaning that a write need not be in constant time anymore.
To ensure the implementation also works when n and/or m is larger than the number of threads d on the GPU, we encapsulate the if -then blocks at lines 1-4, 7-9, 10-15, 16-18, 19-21 and 22-31 of Algorithm 3 each in a for-loop, in which every thread accesses not only the data elements associated with its global index i, but also, if needed, the elements with index i + d, i + 2d, etc., as long as the indices are valid. Table 1 shows the results of the experiments we conducted. The |Blocks| column indicates the number of different blocks at the end of the algorithm, where each block contains only bisimilar states. With #It we refer to the number of while-loop iterations that were executed (see Algorithm 3), before all blocks became stable. The T pre give the preprocessing times in milliseconds, which includes doing the memory transfers to the GPU, sorting the transitions and partitioning. The T alg give the times of the core algorithm, in milliseconds. The T total is the sum of the preprocessing and the algorithm, in milliseconds. We have not included the loading times for the files and the first CUDA API call that initializes the device. We ran each benchmark 10 times and took the averages. The standard deviation of the total times varied between 0% and 3% of the average, thus 10 runs are sufficient. All the times are rounded with respect to the standard error of the mean.
We see that the bound as proven in Section 4.2 (k ≤ 3n) is indeed respected, #It/n is at most 2.20, and most of the time below that. The number of iterations is tightly related to the amount of blocks that the final partition has, the #It/|Blocks| column varies between 1.00 and 2.53. This can be understood by the fact that each iteration either splits one or more blocks or marks a block as stable, and all blocks must be checked on stability at least once. This also means that for certain LTSs the algorithm scales better than linearly in n. The preprocessing often takes the same amount of time (about a few milliseconds). Exceptions are those cases with a large number of actions and/or transitions.
Concerning the GPU run times, it is not true that each iteration takes the same amount of time. A GPU is not a perfect PRAM machine. There are two key differences. Firstly, we suspect that the algorithm is memory bound since it is performing a limited amount of computations. The memory accesses are irregular, i.e., random, which caches can partially compensate, but for sufficiently large n and m, the caches cannot contain all the data. This means that as the LTSs become larger, memory accesses become relatively slower. Secondly, at a certain moment, the maximum number of threads that a GPU can run in parallel is achieved, and adding more threads will mean more run time. These two effects can best be seen in the column headed by T alg /#It, which corresponds to the time per iteration. The values are around 0.02 up to 300, 000 transitions, but for a higher number of states and transitions, the amount of time per iteration increases.

Experimental comparison
We compared our implementation (Par-BCRP) with an implementation of the algorithm by Lee and Rajasekaran (LR) [13] on GPUs, and the optimized GPU implementation by Wijs based on signature-based bisimilarity checking [2], with multi-way splitting (Wms) and with single-way splitting (Wss) [23]. Multi-way splitting indicates that a block is split in multiple blocks at once, which is achieved in signature-based algorithms by computing a signature for each state in every partition refinement iteration, and splitting each block off into sets of states, each containing all the states with the same signature. The signature of a state is derived from the labels of the blocks that this state can reach in the current partition. The running times of the different algorithms can be found in Table 2. Similarly to our previous benchmarks, the algorithms were run 10 times on the same machine using the same VLTS benchmark suite with a time-out of 10 minutes. In some cases, the non-deterministic behaviour of the algorithms Wms and Wss led to high variations in the runs. In cases where the standard error of the mean was more than 5% of the mean value, we have added the standard error in Table 2 in between parentheses. Furthermore, all the results are rounded with respect to the standard error of the mean. As a pre-processing step for the LR, Wms and Wss algorithms the input LTSs need to be sorted. We did not include this in the times, nor the reading of files and the first CUDA API call (which initializes the GPU).
This comparison confirms the expectation that our algorithm in all cases (except one small LTS) out-performs LR. This confirms our expectation that LR is not suitable for massive parallel devices such as GPUs. Furthermore, the comparison teaches that in most cases our algorithm (Par-BCRP) outperforms Wss. In some benchmarks (Cwi 1 2, Cwi 214 684, Cwi 2165 8723 and Cwi 2416 17605) Wss is more than twice as fast, but in 16 other cases our algorithm is more than twice as fast. The last comparison shows us that our algorithm does not out-perform Wms. Wms employs multi-way splitting which is known to be very effective in practice. Contrary to our implementation, Wms is highly optimized for GPUs while the focus of the current work is to improve the theoretical bounds and describe a general algorithm.
In order to understand the difference between Wms and our algorithm better, we analysed the complexity of Wms [23]. In general this algorithm is quadratic in time, and the linearity claim in [23] depends on the assumption that the fan-out of 'practical' transition systems is bounded, i.e., every state has no more than c outgoing transitions for c a (low) constant.
We designed the transition systems Fan out n for n ∈ N + to illustrate the difference. The This LTS has n states, 3n − 3 transitions and a maximum out degree of n transitions. Figure 3 shows the results of calculating the bisimulation equivalence classes for Fan out n , with Wms and Par-BCRP. It is clear that the run time for Wms increases quadratically as the number of states grows linearly, already becoming untenable for a small amount of states. On the other hand, in conformance with our analysis, our algorithm scales linearly.

Weaker PRAM models
Algorithm 2 relies on concurrent writes to perform the constant time leader election and the choice of splitter. This means that the algorithm does not work on a weaker PRAM model. In this section we describe a modification for the common CRCW PRAM and a limitation for the ERCW PRAM.
It is shown in [12] that any priority CRCW PRAM using n processors and m memory cells can be simulated by a common CRCW PRAM with O(n 2 ) processors and O(m 2 ) memory cells. For our problem, a common CRCW PRAM with O(n 2 ) processors and no extra memory can solve leader election.
This leader election on the common CRCW PRAM is given in Algorithm 6. Every processor is indexed as P i,j for all i, j ∈ {0, . . . , n − 1} for exactly n 2 processors. First, if P i,j has a state with index i that is eligible to be the leader of a new block (line 1), it writes block i , i.e., the index of the block the state is currently a member of, to position i in a list new leader . In the next step, P i,j replaces new leader j with 0 if new leader i = new leader j and i < j. In other words, if P i,j encounters two states that can become the new leader, it selects the one with the smallest index. This is possibly a concurrent write, but all writes involve the same value 0, hence this is allowed by the common CRCW PRAM. Next, if for P i,j , new leader i = 0, it writes the value i to new leader block i at line 7. For a given block block i , the condition at line 6 only holds for the state with the largest index among the states that are split from block i , hence there is at most one value is written.
Leader election on the ERCW PRAM is not possible in constant time, which follows from a result by Cook et al. [6,Thm 4.]. This result says that all functions that have a critical input are in Ω(log n) on ERCW PRAMs. A bit sequence I of size n is critical for a function f : {0, 1} n → {0, 1} iff for any I obtained by flipping exactly one bit in I we have f (I) = f (I ). Leader election can be seen as a function f : S → {0, 1}, where f (i) = 1 iff i is elected as a new leader. This function has a critical input, namely the to be chosen leader.

Related work
In [13] Lee and Rajasekaran study RCPP. They implement a parallel version of Kanellakis-Smolka that runs in O(n log n) time on m log n log log n CRCW PRAM processors. In [17] they present a different algorithm based on Paige and Tarjan's algorithm [16] that has the same run time of O(n log n) but using only m n log n CREW processors. Jeong et al. [10] presented a linear time parallel algorithm, but it is probablistic in the sense that it has a non-zero chance to output the wrong result. Furthermore, Wijs [23] presented a GPU implementation of an algorithm to solve the strong and branching bisimulation partition refinement problems but although efficient for many practical cases, it has a quadratic time complexity.
In a distributed setting, Blom and Orzan studied algorithms for refinement [2]. Those algorithms use message passing as ways of communication between different workers in a network and rely on a small number of processors. Therefore, they are very different in nature than our algorithm. Those algorithms were extended and optimized for branching bisimulation [3].

Conclusion
We proposed and implemented an algorithm for RCPP and BCPP. We proved that the algorithm stops in O(n + |Act|) steps on max(n, m) CRCW PRAM processors. We implemented the algorithm for BCPP in CUDA, and conducted experiments that show the potential to compute bisimulation in practice in linear time. Further advances in parallel hardware will make this more feasible.
For future work, it is interesting to investigate whether RCPP can be solved in sublinear time, that is O(n ) for a < 1, as requested in [13]. It is also intriguing whether the practical effectiveness of the algorithm in [23] by splitting blocks simultaneously can be combined with our algorithm, while preserving the linear time upperbound. Furthermore, it remains an open question whether our algorithm can be generalised for weaker bisimulations, such as weak and branching bisimulation [22,9]. The main challenge here is that the transitive closure of so-called internal steps needs to be taken into account.