IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

1

Parallelizing Exact and Approximate String
Matching via Inclusive Scan on a GPU
Yasuaki Mitani, Fumihiko Ino, Member, IEEE, and Kenichi Hagihara
Abstract—In this study, to substantially improve the runtimes of exact and approximate string matching algorithms, we propose atribrid parallel method for bit-parallel algorithms such as the Shift-Or and Wu-Manber algorithms. Our underlying idea is to interpret
bit-parallel algorithms as inclusive-scan operations, which allow these bit-parallel algorithms to run efficiently on a graphics processing
unit (GPU); we achieve this speed-up here because inclusive-scan operations not only eliminate duplicate searches between threads
but also realize a GPU-friendly memory access pattern that maximizes memory read/write throughput. To realize our ideas, we first
define two binary operators and then present a proof regarding the associativity of these operators, which is necessary for theparallelization of the inclusive-scan operations. Finally, we integrate the inclusive-scan scheme into a previous segmentation-based
scheme to maximize search throughput, identifying the best tradeoff point between synchronization cost and duplicate work. Throughour experiments, we compared our proposed method with previous segmentation-based methods and indexing-based sequencealigners. For online string matching, our proposed method performed 6.7–16.7 times faster than previous methods, achieving a searchthroughput of up to 1.88 terabits per second (Tbps) on a GeForce GTX TITAN X GPU. We therefore conclude that our proposed
method is quite effective for decreasing the runtimes of online string matching of short patterns.
Index Terms—String matching, bit-parallel algorithm, inclusive scan, Shift-Or algorithm, Wu-Manber algorithm, GPU.

F

1I NTRODUCTION

T

HE problem of string matching is to find all locations
at which pattern P of length m matches a substring of
text, with T being of length n. Exact string matching finds
all occurrences, if any, of pattern P in text T . Conversely,
approximate string matching finds all strings whose edit
distance from P is smaller than error k , i.e., the tolerable
maximum edit distance. String matching is a fundamental
problem that occurs in a wide range of practical applications, including computational biology [1], information
retrieval [2], and network intrusion detection [3].
String matching algorithms have been intensively studied for over 40 years using several approaches, including
dynamic programming, automata, bit parallelism, filtering,and indexing approaches [4]. Among these, automata-basedalgorithms have the best worst-case time, i.e., O(n) or linear
time, which is therefore the lower bound of this irregularproblem [4]. These automata-based algorithms can be further classified into two groups according to their underlying
automata, namely, deterministic finite automata (DFAs) or
nondeterministic finite automata (NFAs). For DFAs, theKnuth-Morris-Pratt [5] and Aho-Corasick [6] algorithms
have been widely used for exact string matching problems.
With respect to NFAs, a variety of bit-parallel approaches
have been used in practice as fast algorithms for exact andapproximate string matching problems [7]–[9]. In [7], [8],Baeza-Yates and Gonnet proposed the Shift-Or (SO) algorithm, which finds exact matching positions in O(nm/w)••

Y. Mitani is with the Development Head Office, DWANGO Corporation,Ltd., 4-12-15 Ginza, Chuo-ku, Tokyo 104-0061, Japan.
F. Ino and K. Hagihara are with the Graduate School of Information
Science and Technology, Osaka University, 1-5 Yamadaoka, Suita, Osaka
565-0871, Japan. E-mail: ino@ist.osaka-u.ac.jp

Manuscript received on 23 June, 2016; revised 27 Nov. 2016.

time, where w is the size of a computer word. Although
this complexity is not optimal, the SO algorithm is usually
fast in practice, particularly for relatively short patterns
that fit into a computer word (i.e., m ≤ w); this inherent
speed exists because the algorithm takes advantage of the
intrinsic parallelism of bit operations inside a computer
word. Exploiting this bit parallelism reduces the number ofoperations by a factor of up to w. The SO algorithm also
has advantages in terms of data locality because its simple
data structure regularizes memory access patterns. In [9],Wu and Manber extended this SO algorithm to solve theapproximate string matching problem in O(nm/wk) time.
In addition to bit parallelism, other data-parallel methods have been presented to achieve further acceleration
on multicore CPUs [10], graphics processing units (GPUs)[11], and Xeon Phi coprocessors [12]. These hybrid methods exploit both bit parallelism and data parallelism bypartitioning the text into smaller segments, which are thenprocessed in parallel by running a bit-parallel algorithm
for each segment. For example, GBPR [11] deployed thishybrid approach to realize a high search throughput ofapproximately 75 gigabits per second (Gbps) on a GeForceGTX 680 GPU, which was seven to 11 times faster thanthat of a multithreaded CPU implementation. Here searchthroughput ρ is given by ρ = n/t, where t is the GPU
execution time spent for the search.
Nonetheless, using a segmentation-based scheme canresult in a degradation of efficiency on a massively parallelmachine because halo regions of size m − 1 are required tocorrectly find a match ending at the head of a segment,
as depicted in Fig. 1a. These halo regions add substantial
overhead because the characters within the halo region must
be read twice, i.e., once for each of the adjacent segments.
Given this duplicate searching, using a segmentation-based

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X
n

n
Text

Text

s

Segmentation

Text

Parallel prefix-sum
Data-parallel search

Halo region m 1
Text
Thread 1

Parallel scan
Text

+

+
+

+

+
+
+

Thread 2
Thread 3(a)Thread 4

Thread 4(b)Fig. 1: Parallel string search using (a) a segmentation-based
scheme and (b) an inclusive-scan scheme. The former processes a sequential search algorithm for each segment ofsize s; thus halo regions of size m − 1 are required tofind a potential match ending prematurely at the head ofa segment. In contrast, the latter eliminates halo regions but
requires global synchronization to obtain the entire scan.
scheme results in lower efficiency with strong scaling, where
an increasing number of processing elements (i.e., segments)
are applied to a fixed-size problem (i.e., text). Furthermore,parallel threads suffer from strided memory accesses when
they simultaneously access different segments. Thus, thechallenge of realizing a highly efficient solution remains
for increasing the runtimes of bit-parallel algorithms on a
massively parallel machine; more specifically, we need alow-cost emerging accelerator capable of exploiting finegrained parallelism with millions of threads.
A prospective approach for eliminating halo regions
is to interpret bit-parallel algorithms with scan primitives[13], which have been used to parallelize many key operations that seem inherently sequential [14], including
sorting, stream compaction, and histogram computation.
With respect to string comparison, Khajeh-Saeed et al. [15]presented an inclusive-scan-based approach for the SmithWaterman algorithm [16], which can be regarded as a general form of the Wu-Manber (WM) algorithm [9] mentioned
above; the Smith-Waterman algorithm computes the similarity scores of arbitrary regions of two strings according
to a flexible scoring system (i.e., substitution matrix andgap penalties). Thus, scan operations can be efficiently computed on the GPU [14], [17] because this type of computation has high locality and regularity in terms of memoryreferences. However, as far as we know, scan primitives
have not been applied to bit-parallel algorithms, whichhave different computing strategies and data structures ascompared to the Smith-Waterman algorithm.
In this paper, we therefore propose a tribrid parallel
method that takes advantages of scan-based parallelization, segmentation-based parallelization, and bit-level parallelization. The scan scheme allows massively parallel
threads to run the SO algorithm [7], [8] to find matchingpatterns without halo regions, as depicted in Fig. 1b; more
specifically, a string matching problem can be interpreted as
an inclusive-scan problem according to Blelloch’s reduction[18]. To enable this interpretation, we find a companion
operator [18] to rewrite the SO algorithm with a recurrence
relation based on associative operators. This associativity
allows the SO algorithm to be parallelized in the samemanner as that of the parallel scan algorithm [19]. Thus, the2automata can be correctly updated by reading charactersat arbitrary positions of the text, meaning that there is noneed to serially read characters from the head of a text
segment. Because the inclusive-scan scheme requires global
synchronization of threads, this scheme is cascaded with asegmentation-based scheme, which avoids synchronization,
but requires halo regions. We also show that scan-based
parallelization can be applied to the WM algorithm forapproximate string matching. We implemented our proposed methods on a GPU compatible with the computeunified device architecture (CUDA) [20]. The developed
implementation, which we called cuShiftOr, is available at
http://www-hagi.ist.osaka-u.ac.jp/research/code/.
In addition to this introduction, the remainder of thispaper is structured as follows. In Section 2, we introduce
related studies regarding approaches to improving the runtimes of various string matching algorithms. In Section 3, we
present an overview of bit-parallel algorithms, which forms
the basis of our proposed method. In Section 4, we describe
our proposed methods for exact and approximate string
matching. Next, in Section 5, we present key implementation issues that must be solved to achieve CUDA-based
acceleration. In Section 6, we show our experimental results
obtained using the Maxwell GPU [21]. Finally, in Section 7,
we conclude our paper and discuss avenues of future work.

2

R ELATED W ORK

Several parallel implementations for accelerating bitparallel algorithms on a GPU have been presented, someof which are detailed below. Previous methods typicallyparallelize a string search algorithm via a segmentationbased scheme. By contrast, cuShiftOr focuses on the associativity of search operations, which enables text to beread in an arbitrary order without halo regions. As such,
cuShiftOr achieves an efficient memory access pattern forthis irregular problem, yielding Tbps speeds for in-coresearch throughput of short patterns of up to 64 characters.
Lin et al. presented GBPR [11], which deploys asegmentation-based scheme to exploit data parallelism of
the WM algorithm with optimization to the hierarchicalmemory architecture of the GPU. Given segment size s,
text of length n is divided into ⌈n/s⌉ overlapping segments
to allow ⌈n/s⌉ threads to independently run the sequential
WM algorithm for their responsible segments. Furthermore,
GBPR reduces the amount of off-chip memory access by utilizing on-chip shared memory [20] as a software-managed
cache. The achieved search throughput of 75 Gbps was
2.8–104.8 times that of various state-of-the-art approaches[22]–[25]. However, this throughput can be increased if sstrided accesses from a warp [20] can be coalesced into asingle access to a contiguous memory region. Here, a warp
is a collection of 32 GPU threads that execute the same
instruction at every clock cycle.
In [26], Tran et al. presented XBitPar, which allows a
GPU to process the WM algorithm with a long pattern of
up to 1024 characters. Their basic idea was to simulate along computer word by letting each thread within a warp
take a 32-bit word as an operand. This simulation is usefulfor dealing with long patterns, such as biological sequences,which usually cannot fit into a 32-bit computer word. By

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

contrast, our proposed method, which deploys a 32/64bit word, can be faster if a 1024-bit word is wasteful forthe input pattern. This waste issue is not critical for multipattern matching problems because the automata states
for multiple patterns can easily be packed into a computer
word. Further, in [10], Kusudo et al. accelerated the SO
algorithm for multipattern matching on a multicore CPU.
Their implementation utilized advanced vector extensions(AVX) [27], which is an instruction set for single-instruction,
multiple-data (SIMD) processing. In [24] and [28], similar
approaches were implemented on a GPU. Given these previous studies, we conclude that scan-based parallelization canbe applied to the family of bit-parallel algorithms designed
for multipattern matching problems.
Similar to the WM algorithm, Baeza-Yates et al. [29]realized approximate string matching by packing automaton states in a different manner. The difference here wasthat they packed the automaton states along the diagonals
instead of along the rows or columns, thereby achieving
O(n) search for short patterns, i.e., mk = O(w). Theiralgorithm was implemented on an NVIDIA Tesla S1070
GPU by Onsjö et al. [25], who achieved a search throughput
of 158 Mbps. Similar to XBitPar [26], they realized a solutionthat supported a long computer word via the warp-based
search technique described above.
Dynamic programming has also been applied to string
search problems. Using dynamic programming, a matrix
called the dynamic programming matrix must be filled
to find matching positions, with new values within thematrix building upon previous values of the matrix. In [30],Myers efficiently parallelized the computation of the dynamic programming matrix by considering the differences
along columns instead of the values within the columns.
He exploited bit parallelism in the matrix computation,achieving a runtime of O(nm/w) for arbitrary k . Given
its time complexity is independent of k here, the Myers
algorithm is suitable for approximate string matching,In [31], Langner parallelized the Myers algorithm usinga GPU by simultaneously processing matrix elements on thesame antidiagonal. Although this parallelization scheme isefficient for large n and m, efficiency degrades for short patterns because the degree of parallelism is limited to at mostmin(n, m/w) elements located on the same antidiagonal. A
similar GPU-based dynamic programming implementation
was presented by Man et al. [32], who achieved a searchthroughput of 80 Mbps on an NVIDIA GeForce GTX 580.
Lin et al. [33] parallelized the Aho-Corasick algorithm[6] on a GPU. Their implementation, called the parallel
failureless Aho-Corasick (PFAC) algorithm, modifies the
DFA to reduce the number of necessary warp divergences[20], which significantly decreases the efficiency of SIMD
execution. A search throughput of up to 143 Gbps wasachieved on a GeForce GTX 580 for 1998 exact patterns.
This throughput was 14.7 times faster than that of an
OpenMP-based CPU implementation; however, the PFAC
algorithm suffers from load imbalance issues if too many
partial matches are found during searching, as is evident
by throughputs ranging from 11.0–218.7 Gbps. By contrast,
bit-parallel algorithms avoid such load imbalance issues,
opening the door to stable search throughputs regardlessof the number of partial matches.

3

A

T

A

G
0 error

0,0

A
0

T
1

A
2

G
3

2 errors

(a)

A

T

A

G

A

T

A

G

4

1 error2 errors(b)

Fig. 2: NFA examples used for bit-parallel algorithms. Given
pattern ATAG, the SO algorithm [7], [8] generates (a) an
NFA to be used for exact string matching, while (b) theWM algorithm [9] extends the NFA to perform approximate
string matching with at most k errors (in this example,
k = 2). Note that unlabeled transitions match any arbitrary
character.
Finally, an important application category of string
matching is biological sequence alignment. Many exact andapproximate algorithms [34]–[38] have been proposed withindexing architectures. For short patterns, indexing-based
approaches are usually computationally efficient versus indexless (i.e., brute-force) approaches; however, indexingbased approaches assume offline searching and require
extra memory space to maintain indexes. We compare
cuShiftOr with state-of-the-art aligners [34]–[38] later in
Section 6.4.

3

P RELIMINARIES : B IT-PARALLEL A LGORITHMS

In the following subsections, we use a C-like notation forbitwise operators, i.e., let ≪, &, |, and ∼ be the left shift,AND, OR, and NOT operators, respectively. Further, let Z
and F be the sets of non-negative integers and a binary
finite field, respectively. Also, let Fw be the set of binary
bit-vectors of size w.
3.1

Shift-Or Algorithm for Exact Matching

In this subsection, we describe the SO algorithm for exact string matching. Let Σ be an alphabet of characters.
Let T = t1 t2 · · · tn be the text of length n ∈ Z and let
P = p1 p2 · · · pm be the pattern of length m ∈ Z, whereti ∈ Σ (1 ≤ i ≤ n) and pj ∈ Σ (1 ≤ j ≤ m) are the i-thand j -th characters of the text and pattern, respectively. The
SO algorithm [7], [8] accelerates string matching by runningm string comparators in parallel, each reading the same text
position concurrently but with a different starting position.
Each comparator simulates an NFA to express a matchingstate; the current state of the NFA is updated by readinga character from the text. Figure 2a shows an example of
the NFA generated for pattern ATAG. As shown in the
figure, state ri indicates that the first i characters of thepattern exactly match the last i characters of the text, where1 ≤ i ≤ m. Consequently, an exact match ending at tj can bedetected if the NFA reaches accepting state rm after reading
tj , where 1 ≤ j ≤ n.
Bit parallelism can be exploited here by packing thecurrent states of m NFAs into bit-vector R ∈ Fm of size msuch that R contains information on all matches of prefixesof pattern P . Let Rj ∈ Fm be the value of bit-vector
R after tj has been read. Rj is defined as Rj [i] = 0 if

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

p1 p2 · · · pi = tj−i+1 tj−i+2 · · · tj ; otherwise, Rj [i] = 1. Bitvector Rj can then be computed according to a first-orderrecurrence relation{
1m ,if j = 0,(1)
Rj =(Rj−1 ≪ 1) | Mj , if j < 0 ≤ n,
where 1m is an all-one vector of size m and Mj ∈ Fm is a
bit-vector that depends on tj : for all 1 ≤ i ≤ m, Mj [i] = 0,
if tj = pi ; otherwise, Mj [i] = 1. In other words, Mj is amask that can be stored in a table indexed by alphabet Σ.
The SO algorithm can be rapidly processed if a bit-vector
fits within a computer word (i.e., m ≤ w).
3.2

Wu-Manber Algorithm for Approximate Matching

The basic idea behind the WM algorithm for approximate
string matching is to extend the NFA of the SO algorithm
to allow transitions between states associated with error
distance d (0 ≤ d ≤ k). Let rd,i be the state in which thefirst i characters of the pattern match the last i charactersof the text with d errors. In this way, the extended NFA
considers all possible matches with up to k errors. Figure 2b
illustrates an NFA example for pattern ATAG with k = 2.
Similar to the SO algorithm, the WM algorithm [9]exploits bit parallelism by packing the current states of m
NFAs, i.e., bit-vector R ∈ Fm of size m contains informationon all matches of prefixes of pattern P with up to k errors.
Let Rd,j ∈ Fm be the value of the bit-vector after tj has
been read. For approximate matching, we must considertransitions not only with matches but also from insertions,deletions, or substitutions, which we express as


1m ,





(Rd,j−1 ≪ 1) | Mj ,



m−d d


0 ,
1
Rd,j = ((Rd,j−1 ≪ 1) | Mj )
& Rd−1,j−1





&
(Rd−1,j−1 ≪ 1)
 & (Rd−1,j ≪ 1),

3.3

if d = 0, j = 0,
if d = 0, 0 < j ≤ n,
if 0 < d ≤ k, j = 0,

if 0 < d ≤ k, 0 < j ≤ n.
(2)

Properties of Fundamental Operators

To facilitate the understanding of our proposed method, we
describe properties of fundamental bitwise operations via
the lemmas below.

4

4

P ROPOSED PARALLEL M ETHOD

The fundamental goal of our proposed method is to introduce associative binary operator • to translate Eqs. (1) and(2) into recurrence relation{
a0 ,if j = 0,
Rj =(3)
Rj−1 • aj , if j > 0.
To achieve this, we find a companion operator [18] for Eq. (1)and show a proof based on Blelloch’s reduction [18]; byusing the companion operator, a first-order recurrence witha semiassociative operator [18] and an associative operator
can be reduced into Eq. (3). We also extend this idea for
Eq. (2), i.e., a recurrence with three different operators,≪, |, and &, by defining an associative operator and asemiassociative operator required for Blelloch’s reduction.
4.1

Exact Matching

For exact string matching problems, we introduce bitwise
operator † as the key operator •.
Definition 1. Bitwise operator † : (Z × Fw ) × (Z × Fw ) →Z × Fw for pairs of a non-negative integer and a bit-vector of size
w is defined asdef⟨u1 , x1 ⟩ † ⟨u2 , x2 ⟩ = ⟨u1 +u2 , (x1 ≪ u2 ) | x2 ⟩,(4)

where u1 , u2 ∈ Z and x1 , x2 ∈ Fw . The operator + here is thecompanion operator of the semiassociative operator ≪.
Using associative operator †, an exact string matchingproblem can be solved as an inclusive-scan problem. Algorithm 1 shows our proposed method for exact matching,∑jwhere we use
l=0 Al to denote a summation over †:∑jm
l=0 Al = A0 †A1 †· · ·†Aj , where A0 , A1 , . . . , Al ∈ Z × F .
The basic idea behind operator † is that it pairs up a maskwith an integer to correctly compute accumulated mask
Rj ∈ Z × Fm , thus accumulated mask Rj from Algorithm 1
considers all possible states that can be reached after Tj has
been read. Consequently, the bit-vector can be updated byan inclusive-scan scheme as presented in lines 4–6 of Algorithm 1. Note that this algorithm assumes m ≤ 64 (= w) toachieve acceleration on a GPU.
Algorithm 1 SO exact matching with parallel inclusive-scan
Input: Text T of n characters and pattern P of m characters
Output: Bit sequence X of size n that indicates matchingpositions
1: Initialize the masks M1 , M2 , · · · , Mn according to T

and P

Lemma 3. ≪ is right-distributive over |.

2: A0 ← ⟨0, 1m ⟩3: Aj ← ⟨1, Mj ⟩, for all 1 ≤ j ≤ n
4: for j ← 0 to n do in parallel ▷ Parallel inclusive-scan∑
5:
Rj ← jl=0 Al▷ Summation over †6: end for
7: ⟨uj , xj ⟩ ← Rj , for all 0 ≤ j ≤ n
8: X[j] ← xj [m], for all 1 ≤ j ≤ n
9: return X

The proofs of the above lemmas are presented
in the supplementary material, which can befound on the Computer Society Digital Library at
http://doi.ieeecomputersociety.org/10.1109/TPDS.xxxx.xx.

Next, we show a proof of the correctness of Algorithm 1.
First, we show that † is an associative operator required forparallelizing lines 4–6 of Algorithm 1.

Lemma 1. The bitwise left shift operator ≪: Fw × Z → Fw
satisfies (x ≪ u) ≪ v = x ≪ (u + v), where x ∈ Fw is a
bit-vector of size w and u, v ∈ Z are non-negative integers.
Lemma 2. The bitwise OR operator |: Fw × Fw → Fw is
associative.

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

Theorem 1. † is associative.
Proof. The proof here is given by Lemmas 1, 2, and 3. Let
x1 , x2 , x3 ∈ Fw and u1 , u2 , u3 ∈ Z. Then,(⟨u1 , x1 ⟩ † ⟨u2 , x2 ⟩) † ⟨u3 , x3 ⟩= ⟨u1 +u2 , (x1 ≪ u2 ) | x2 ⟩ † ⟨u3 , x3 ⟩= ⟨(u1 +u2 ) + u3 , (((x1 ≪ u2 ) | x2 ) ≪ u3 ) | x3 ⟩= ⟨u1 +u2 +u3 , (((x1 ≪ u2 ) ≪ u3 ) | (x2 ≪ u3 )) | x3 ⟩= ⟨u1 +u2 +u3 , (x1 ≪ (u2 +u3 )) | (x2 ≪ u3 ) | x3 ⟩. (5)Conversely,

5

Definition 4. Bitwise operator ‡ : (Z × Gw ) × (Z × Gw ) →Z × Gw is defined asdef⟨u1 , ⟨x1 , y1 ⟩⟩ ‡ ⟨u2 , ⟨x2 , y2 ⟩⟩ =⟨u1 +u2 , (⟨x1 , y1 ⟩ ≪2 u2 ) ⊎ ⟨x2 , y2 ⟩⟩, (10)where u1 , u2 ∈ Z, and x1 , x2 , y1 , y2 ∈ Fw . The operator + hereis the companion operator of the semiassociative operator ≪2 .
Using operator ‡, the approximate string matching problem can be solved via an inclusive-scan scheme, as shownin Algorithm 2.

⟨u1 , x1 ⟩ † (⟨u2 , x2 ⟩ † ⟨u3 , x3 ⟩)= ⟨u1 , x1 ⟩ † ⟨u2 +u3 , (x2 ≪ u3 ) | x3 ⟩= ⟨u1 + (u2 +u3 ), (x1 ≪ (u2 +u3 )) | ((x2 ≪ u3 ) | x3 )⟩= ⟨u1 +u2 +u3 , (x1 ≪ (u2 +u3 )) | (x2 ≪ u3 ) | x3 ⟩. (6)Algorithm 2 WM approximate matching with parallel
inclusive-scan
Input: Text T of n characters and pattern P of m characters
Output: Bit sequence Y of size n that indicates approximate
matching positions with at most k errors

Thus, Eq. (5) is equivalent to Eq. (6), and thereby † is
associative.

1: Initialize the masks M1 , M2 , · · · , Mn according to T

Next, we show that Algorithm 1 produces the same
results of the SO algorithm.
Theorem 2. Let ⟨u0 , x0 ⟩, ⟨u1 , x1 ⟩, . . . , ⟨un , xn ⟩ ∈ Z × Fm be
the inclusive scans of pairs obtained at line 7 of Algorithm 1,
where m ∈ Z. Then, for all 0 ≤ j ≤ n, xj = Rj , where
Rj ∈ Fm is the j -th bit-vector obtained by Eq. (1).
Proof. According to Algorithm 1, j -th inclusive scan⟨uj , xj ⟩, where 0 ≤ j ≤ n, can be expressed as recurrence
relation{⟨0, 1m ⟩,if j = 0,⟨uj , xj ⟩ =(7)⟨uj , xj−1 ⟩ † ⟨1, Mj ⟩, if 0 < j ≤ n.
Conversely, Eq. (4) gives ⟨uj , xj−1 ⟩ † ⟨1, Mj ⟩ = ⟨uj +1, (xj−1 ≪ 1) | Mj ⟩, meaning that Eq. (1) is part of Eq. (7);the second component of the pairs in Eq. (7) is equal to
Eq. (1); hence, xj = Rj , for all 0 ≤ j ≤ n.
4.2

Approximate Matching

We cannot directly apply our exact matching approach toan approximation problem because Eq. (2) consists of | and&, which corrupt the associative property. As an exampleof this, (1 & 0) | 1 ̸= 1 & (0 | 1). To address this issue,
we introduce three additional bitwise operators ≪2 , ⊎, and‡, where ‡ is the target associative operator • constructedusing ≪2 and ⊎. These operators are defined below.
Definition 2. Bitwise left shift operator ≪2 : Gw × Z → Gw fordef
pairs of bit vectors, where Gw = Fw × Fw , is defined as
def⟨x, y⟩ ≪2 u = ⟨x ≪ u, y ≪ u⟩,(8)

where x, y ∈ Fw and u ∈ Z. The operator ≪2 here is thesemiassociative operator required for Blelloch’s reduction.
Definition 3. Bitwise operator ⊎ : Gw ×Gw → Gw is defined asdef⟨x1 , y1 ⟩ ⊎ ⟨x2 , y2 ⟩ = ⟨x1 | x2 , (y1 & ∼x2 ) | y2 ⟩,(9)where x1 , x2 , y1 , y2 ∈ Fw . The operator ⊎ here is the associative
operator required for Blelloch’s reduction.

and P

2: Compute R0,j , for all 0 ≤ j≤ n using lines 2–6 ofAlgorithm 1
3: ⟨u0,j , y0,j ⟩ ← R0,j , for all 0 ≤ j ≤ n
4: for d ← 1 to k do5:
Nj ← yd−1,j−1 & (yd−1,j−1 ≪ 1) & (yd−1,j ≪ 1),
for all 1 ≤ j ≤ n6:Ad,0 ← ⟨0, ⟨0m , 1m−d 0d ⟩⟩7:Ad,j ← ⟨1, ⟨∼Nj , (Mj & Nj )⟩⟩, for all 1 ≤ j ≤ n
8:for j ← 0 to n do in parallel▷ Parallel
inclusive-scan ∑9:Rd,j ← jl=0 Ad,l
▷ Summation over ‡10:
end for11:⟨ud,j , ⟨xd,j , yd,j ⟩⟩ ← Rd,j , for all 0 ≤ j ≤ n
12: end for
13: Y [j] ← yk,j [m], for all 1 ≤ j ≤ n
14: return Y
Next, we show a proof of the correctness of Algorithm 2.
Equations (1) and (2) indicate that R0,j , where 0 ≤ j ≤ n,can be obtained by processing the SO algorithm, i.e., lines 2
and 3 of Algorithm 2 correctly produce exact matching
results for d = 0. Regarding remaining lines 4–12, which
consider cases of 0 < d ≤ k , we show below that Algorithm 2 produces the same results as the WM algorithm.
First, we describe some of the properties of the newly
introduced operators. Note here that the associativity of ‡can be shown in a similar manner as that of †, because Eq. (4)for † can be mapped to Eq. (10) for ‡ by substituting ≪ and| for ≪2 and ⊎, respectively. Accordingly, we show somefundamental properties of ≪2 and ⊎ below.
Lemma 4. ≪2 has the property(⟨x, y⟩ ≪2 u) ≪2 v = ⟨x, y⟩ ≪2 (u + v),
where x, y ∈ Fw and u, v ∈ Z.
Proof. The proof here can easily be shown in the samemanner as that of Lemma 1.
Lemma 5. ⊎ is associative.
Proof. When w = 1, (⟨x1 , y1 ⟩ ⊎ ⟨x2 , y2 ⟩) ⊎ ⟨x3 , y3 ⟩ =⟨x1 , y1 ⟩⊎(⟨x2 , y2 ⟩⊎⟨x3 , y3 ⟩), for all x1 , x2 , x3 , y1 , y2 , y3 ∈ F.

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

This equivalence can easily be confirmed as presented in the
supplementary material. This associativity also exists whenw > 1 because arbitrary bits of the output vectors can becomputed independently of one another.
Lemma 6. ≪2 is right-distributive over ⊎.
Proof. The proof here can easily be shown in the samemanner as that of Lemma 3.

We next extend Theorem 4 to handle inclusive scansobtained with ≪, i.e., in addition to & and |.
∪Theorem 5. Let a0 , a1 , · · · , an ∈ F Fm be a sequence ofn + 1 bit-vectors of size 1 and m ∈ Z. Consider inclusive scans
R0 , R1 , . . . , Rn of these bit-vectors that are summed up using
Eq. (11) with an extended set of operators, opj ∈ {&, |, ≪}. Let⟨xj , yj ⟩ be a pair of bit-vectors given by recurrence relation⟨xj , yj ⟩ = m⟨0 , a0 ⟩,



⟨x
m
j−1 , yj−1 ⟩ ⊎ ⟨0 , aj ⟩,⟨xj−1 , yj−1 ⟩ ⊎ ⟨∼aj , 0m ⟩,⟨xj−1 , yj−1 ⟩ ≪2 aj ,Second, we show that ‡ is an associative operator, required for parallelizing lines 8–10 of Algorithm 2.
Theorem 3. ‡ is associative.
Proof. The proof here can be shown in the same manner asthat of Theorem 1 by using Lemmas 4, 5, and 6.
Third, we show that arbitrary inclusive-scans obtainedwith &, |, and/or ≪ can be expressed with ⊎ and ≪2
according to a rewrite rule. To facilitate this process, we firstconsider the first two operators, & and |.
Theorem 4. Let a0 , a1 , · · · , an ∈ Fm be a sequence of n bitvectors of size m ∈ Z. Consider inclusive scans of these bit-vectors
summed up with & and/or | as{a0 ,if j = 0,
Rj =(11)
Rj−1 opj aj , if 0 < j ≤ n,
where opj ∈ {&, |} and aj , Rj ∈ Fm . Next, let ⟨xj , yj ⟩, where
0 ≤ j ≤ n, be a pair of bit-vectors given by recurrence relation⟨xj , yj ⟩ =m
if j = 0,⟨0 , a0 ⟩,⟨xj−1 , yj−1 ⟩ ⊎ ⟨0m , aj ⟩,if 0 < j ≤ n, opj =|,m⟨xj−1 , yj−1 ⟩ ⊎ ⟨∼aj , 0 ⟩, if 0 < j ≤ n, opj =&,(12)
where xj , yj ∈ Fm . Then, there is a rewrite rule that maps
Eq. (11) to Eq. (12), i.e., Rj = yj , for all 0 ≤ j ≤ n.
Proof. We show a proof for the case of m = 1 by mathematical induction on n. When n = 1, the statement obviously
holds. Assume that the statement holds for n = k ′ (> 1): forall 0 ≤ j ≤ k ′ , Rj = yj , where Rj , yj ∈ F. Let n = k ′ + 1.
According to Eq. (11),Rk′ , if ak′ = 0, opk′ =|,



0,
if ak′ = 0, opk′ =&,Rk′ +1 =
(13)

1,
if
ak′ = 1, opk′ =|,Rk′ , if ak′ = 1, opk′ =& .
Equation (12) indicates that⟨xk′ , yk′ ⟩,
⟨1, 0⟩,⟨xk′ +1 , yk′ +1 ⟩ =⟨xk′ , 1⟩,⟨xk′ , yk′ ⟩,if ak′ = 0, opk′ =|,if ak′ = 0, opk′ =&,if ak′ = 1, opk′ =|,if ak′ = 1, opk′ =& .
(14)

Hence, Rj = yj , for all 0 ≤ j ≤ k ′ + 1; the statement holds
for m = 1. The statement also holds for m > 1 because
arbitrary bits of the output bit-vectors can be computed
independently of one another.

6

if j = 0,
if 0 < j ≤ n, opj =|,if 0 < j ≤ n, opj =&,if 0 < j ≤ n, opj =≪,(15)

where xj , yj ∈ Fm . Then, Rj = yj , for all 0 ≤ j ≤ n.
Proof. Let l be the number of ≪ operators that appear in
Rn = a0 op1 a1 op2 · · · opn an . The proof is given bymathematical induction on l. When l = 0 (i.e., Rn does nothave ≪), the statement holds from Theorem 4. Assume thatthe statement holds for l = k ′ (> 0) and let l = k ′ + 1. Let
J be the largest j such that opJ =≪. Then, the statement forl = k ′ + 1 can be proved in the following three steps.
•••

For all 1 ≤ j < J , Rj = yj holds because RJ−1
includes k ′ ≪ operators.
For j = J , Rj = yj because (1) opJ =≪ indicates
RJ = RJ−1 ≪ aJ and (2) Eq. (15) gives ⟨xJ , yJ ⟩ =⟨xJ−1 , yJ−1 ⟩ ≪2 aJ = ⟨xJ−1 ≪ aJ , yJ−1 ≪ aJ ⟩.
For all J < j ≤ n, opj ̸=≪ holds; that is, opj ∈ {&, |}holds. Hence, Theorem 4 can be applied to Rn =RJ opJ+1 aJ+1 opJ+2 · · · opn an , which indicates
Rj = yj , for all J < j ≤ n.

Thus, the statement holds for arbitrary 0 ≤ l ≤ n, i.e., Rj =yj holds for all 0 ≤ j ≤ n.
Theorem 6. For all 0 < d ≤ k , Rd,j of Eq. (2) is equivalent toyd,j given by⟨xd,j , yd,j ⟩ =

m
m−d d

0 ⟩,
if j = 0,⟨0 , 1(⟨xd,j−1 , yd,j−1 ⟩ ≪2 1)⊎ ⟨∼Nd,j , Mj & Nd,j ⟩, if 0 < j ≤ n,

(16)

where Nd,j = yd−1,j−1 & (yd−1,j−1 ≪ 1) & (yd−1,j ≪ 1).
Proof. When j = 0, it obviously holds. When 0 < j ≤n, Theorem 5 translates Eq. (2) into a recurrence relation
with ⊎ and ≪2 , i.e., when 0 < d ≤ k and 0 < j ≤ n,
Eq. (2) is a recurrence relation with op1 =≪, op2 =|, and
op3 =&. Therefore, applying Theorem 5 to Eq. (2) produces
recurrence relation⟨xd,j , yd,j ⟩= (⟨xd,j−1 , yd,j−1 ⟩ ≪2 1) ⊎ ⟨0m , Mj ⟩ ⊎ ⟨∼Nd,j , 0m )⟩= (⟨xd,j−1 , yd,j−1 ⟩ ≪2 1) ⊎ ⟨∼Nd,j , Mj & Nd,j ⟩.
Thus, Eq. (2) can be mapped to Eq. (16) using ⊎ and ≪2
instead of | and ≪.

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

7

Finally, we next show that Algorithm 2 produces thesame results as that of the WM algorithm.
Theorem 7. Let ⟨ud,j , ⟨xd,j , yd,j ⟩⟩ be an inclusive-scan obtained at line 11 of Algorithm 2. Then, for all 0 ≤ j ≤ n and
0 ≤ d ≤ k , yd,j = Rd,j , where Rd,j ∈ Fm is the bit-vector
obtained by Eq. (2).
Proof. When d = 0 (i.e., exact matching), Theorem 2 provides the proof. We therefore consider the case in whichd > 0. According to Algorithm 2, ⟨ud,j , ⟨xd,j , yd,j ⟩⟩, where
0 ≤ j ≤ n and 0 < d ≤ k , is given by recurrence relation⟨ud,j , ⟨xd,j , yd,j ⟩⟩ =m m−d d
0 ⟩⟩,if j = 0,⟨0, ⟨0 , 1⟨ud,j−1 , ⟨xd,j−1 , yd,j−1 ⟩⟩‡⟨1, ⟨∼Md,j , (Mj & Nd,j )⟩⟩, if 0 < j ≤ n.
Definition 4 rewrites the second components of this equation into Eq. (16), meaning that yd,j = Rd,j for all 0 < d ≤ k
and 0 ≤ j ≤ n.
4.3

Reducing the Time Complexity

The time complexity of the parallel inclusive-scan algorithm[19] is O(n log n), where n is the number of elements to besummed. Therefore, our proposed parallel method incurs aparallelization overhead versus the original SO algorithm,
which finds matching positions in O(nm/w). To reduce thisoverhead, we adopt an optimization strategy that reduces
the time complexity to O(nm/w log m).
We consider here the SO algorithm, but the strategy presented below can easily be applied to the WM algorithm byinterpreting the associative operator † as ‡ for approximate
string matching.
According to Eq. (1), we have

Rj = Mj | (Mj−1 ≪ 1) | · · · | (Mj−m+1 ≪ (m − 1))| (Mj−m+2 ≪ (m − 2)) | · · · | (M0 ≪ j),where 0 < j ≤ n. Because x ≪ k = 0m holds for all x ∈ Fm
and k ≥ m, Eq. (1) can be rewritten as

Rj = Mj | (Mj−1 ≪ 1) | . . . | (Mj−m+1 ≪ (m−1)). (17)This equation specifies that exact inclusive-scans are not required to obtain results of the SO algorithm; in other words,
as the original SO algorithm implies, only the latest m − 1characters in the text are required for detecting an exact
string match ending at the current position. Consequently,∑jline 5 of Algorithm 1, i.e., Rj ← l=0 Al , can be replaced∑jwith Rj ←l=j−m+1 Al to reduce the time complexity
without loss of generality. In this case, n elements are
summed via log m steps, and thereby the time complexity is
reduced from O(n log n) to O(n log m) if m = O(w).
4.4

Reducing Global Synchronization

A major drawback of the parallel scan algorithm [19] is thatit synchronizes all threads, requiring each to walk throughlog m steps. Conversely, the CUDA programming modelprohibits global synchronization within a kernel [20] (i.e.,
GPU program); in other words, a kernel completion implicitly synchronizes all GPU threads. Therefore, the kernel

s

s
m ✄1

… A

s ( i ✆ 2 ) ☎1

Phase 1+…
…

As (i ✠ 2 )✟ s As (i ✞1)✝1

+

+

…
…

+
Phase 2…

Ci ✁1,1

…
…
Ci ✂1

…

+
+

+ +
Ci 1, sAs ( i ✆1)☎ s

Ci ,1

… + +
… C

i,s

Ci

…

Fig. 3: Overview of the disjoint version of our tribrid
method. In the first phase, we perform the parallel inclusivescan algorithm within each disjoint segment of s (≥ m)elements. Next, each segment Ci is updated by summing up
the last element Ci−1,s of its left neighbor Ci−1 , thus taking
m−1 preceding elements into consideration. By contrast, theoverlapping version can be processed within a single phase,in which parallel inclusive-scans for contiguous s + m − 1elements are processed independently.

must be invoked log m times to run Algorithm 1 on theGPU, thus search throughput suffers from the increased
amount of off-chip memory access. Data in register files
and shared memory have lifetimes of the corresponding
thread block [20], thus threads must fetch data from off-chip
memory at the beginning of every kernel invocation.
To reduce this overhead, we adopt a tribrid strategy that
integrates the inclusive-scan scheme into a segmentationbased scheme. The following two segmented variations
adapt the tribrid strategy to CUDA: (1) an overlapping
version that requires halo regions and (2) a disjoint version that eliminates halo regions but requires one global
synchronization instead of log m global synchronizations.
Both versions reduce synchronization costs by replacing
global synchronization with local synchronization. The goal
here is to limit synchronization to within a small segment.
Similar to the segmentation-based scheme illustrated in Fig.
1a, the overlapping version allows different segments tobe processed independently by running the parallel scanalgorithm for s+m−1 contiguous elements. By contrast, thedisjoint version runs the parallel scan algorithm for each disjoint segment of s elements but requires an additional phase
to correctly handle the segmentation border, as illustrated
in Fig. 3. Note that Theorem 2 realizes this disjoint version;
thus, there is a tradeoff between the synchronization costand halo region size (i.e., duplicate searching).
Algorithm 3 shows the disjoint version of our proposed
tribrid method, namely SO exact string matching with disjoint segmented parallel inclusive-scan. First, our proposed
method assumes that the text is partitioned into ⌈n/s⌉
segments, where s (≥ m) indicates segment size. Second,
our method independently executes the parallel inclusivescan algorithm for each segment. Thus, the i-th segment,
where 1 ≤ i ≤ ⌈n/s⌉, produces vector Ci ∈ (Z × Fm )s of
def
local inclusive-scans such that Ci = (Ci,1 , Ci,2 , · · · , Ci,s ),IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

Algorithm 3 SO exact matching with disjoint segmented
parallel inclusive-scan
Input: Text T of n characters, pattern P of m characters,
and segment size s (≥ m)
Output: Bit sequence X of size n that indicates matchingpositions
1: Initialize the masks M1 , M2 , · · · , Mn according to T

and P

2: A0 ← ⟨0, 1m ⟩3: Aj ← ⟨1, Mj ⟩, for all 1 ≤ j ≤ n
4: for i ← 1 to ⌈n/s⌉ do in parallel5:for j ← 1 to∑s do in parallel
6:
Ci,j ← jl=1 As(i−1)+l7:
end for
8: end for
9: C0,s ← A0
10: for i ← 1 to ⌈n/s⌉ do in parallel11:
for j ← 1 to s do in parallel
12:
Rs(i−1)+j ← Ci−1,s † Ci,j13:
end for
14: end for
15: ⟨uj , xj ⟩ ← Rj , for all 1 ≤ j ≤ n
16: X[j] ← xj [m], for all 1 ≤ j ≤ n
17: return X▷ for each segment▷ for each element▷ Summation over †▷ for each segment▷ for each element

where

Ci,j =j∑

As(i−1)+l .

(18)

l=1

The above equation indicates that the first m − 1 elements of
Ci do not include the m − 1 left neighbors. Therefore, theseelements must refer to left segment Ci−1 to obtain the final
results, as shown in Fig. 3. The final results can be obtainedby simply adding last element Ci−1,s of the left segment
instead of m − 1 elements (i.e., line 12 of Algorithm 3)
because Ci−1,s has already summed its left neighbors atline 6. Note that we avoid a conditional branch by letting all
elements participate in this addition because summation of
more than m − 1 preceding elements yields the same result
as that produced from exactly m−1 preceding elements, i.e.,∑j∑j
for all 1 ≤ j ≤ n, l=0 Al = l=j−m+1 Al , as presented in
Section 4.3.

5

CU S HIFTO R : GPU- BASED I MPLEMENTATIONIn this section, we present implementation issues and solutions for achieving full optimization of our proposed
method on a CUDA-compatible GPU. Solutions are presented for the SO algorithm, but they can easily be applied
to the WM algorithm. We consider the Maxwell architecture[21] as the target GPU; accordingly, the word size w is set to
32 or 64 according to pattern size m.
Table 1 shows an overview of parallel schemes deployedin our cuShiftOr implementation. Our tribrid method canbe efficiently implemented with CUDA [20], which shares
instructions among threads of the same warp. Threads
in a warp are naturally synchronized, and thereby, theinclusive-scan scheme runs efficiently at the warp level
without additional synchronization costs. Further, a shuffle

8

TABLE 1: Parallelization schemes deployed for each level ofthe thread hierarchy. Fine-grained parallelism is exploitedby the parallel scan scheme whereas coarse-grained parallelism is exploited by the segmentation-based scheme.
Level

Parallel scheme

Halo Synchronization

Grid
Segmentation-based yes
Thread block Inclusive-scannoWarpInclusive-scannoThreadSequential bit-parallel nono
yes w/ __syncthreads()
yes w/o additional cost
n/a

instruction [20] is useful to exchange data for computing inclusive scans within a warp [39]. The inclusive-scan
scheme is integrated into the disjoint segmented scheme
to avoid duplicate searches at the lower thread level. As
for the thread block level, cuShiftOr deploys the same disjoint segmented scheme, because warps in the same thread
block can be synchronized without a kernel completion anddata exchange within a thread block is available throughon-chip shared memory [20], namely, a software-managed
cache. By contrast, the overlapping scheme (i.e., with halo
regions) independently runs at the upper level of the thread
hierarchy—i.e., the grid level—to avoid data exchange between different thread blocks. The search process then completes with a single kernel invocation.
In summary, a thread block is responsible for a segment and its corresponding halo region. This responsiblearea is further decomposed into disjoint subsegments, witheach subsegment assigned to a warp via the disjoint segmented scheme described above. Because a warp contains
32 threads, the segment/subsegment size must be a multipleof 32 to maximize resource utilization. To realize this, thehalo size is set to the maximum word size (i.e., 64) for
arbitrary pattern lengths.
As for the data structure, cuShiftOr assumes that a character is stored as a 1-byte unsigned integer; thus, alphabetΣ is allowed to have 256 symbols. Given n-byte text data
and m-byte pattern data, our proposed implementation
produces n-byte data containing bit sequence X ; each bitof X is stored as 1-byte data to simplify the memory
access pattern. These input/output data are stored in offchip global memory [20], where memory read/write (R/W)
throughput can be maximized when each warp accesses
a contiguous memory region of 128 bytes [20]. For suchcontiguous access, memory transactions from threads ofthe same warp can be coalesced into a single transaction.
Further, shared memory holds a table indexed by alphabet
Σ. This table is useful for computing the mask Mj , where
1 ≤ j ≤ n, according to input character tj , as described in
Section 3.1. Finally, we store current bit-vector Rj , where
0 ≤ j ≤ n, in a register to take advantage of intrinsic
parallelism of bit operations.
5.1

Algorithm Cascading for Coalesced AccessOur cuShiftOr implementation realizes memory access coalescing for 1-byte data by assigning α (≥ 4) contiguous characters of text to a thread. Thus, α represents the task granularity for threads. To enable this assignment, the sequential
bit-parallel algorithm is cascaded to the parallel inclusivescan algorithm, as shown in Fig. 4. To more easily describe

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

… A ✁ ✄ ✂✁As✝ (i ✆ 2 ) ☎1

Phase 1.
Local
prefix-sum(per-thread)As✪ (i ✬ 2 ) ✫ 2✪

s (i 2 )++…

As✎(i ✍ 2 )✌1

Phase 2.
Partialprefix-sum(per-warp)…

As✕(i ✔1) ✓ s+sPhase 1 and 2 of Fig. 3
Ci✜✛1,1Phase 3.
Local
addition(per-thread)

…

+

s

+

…

Ci ✥1,1…

+++

Ci ✗1,✖

Ci ✚1, 2✙

Fig. 4: Overview of algorithm cascading. First, each thread
sequentially computes local inclusive-scans of α consecutive elements. Second, inclusive scans for each α value arecomputed by a warp in parallel according to Algorithm 3.
Finally, we obtain global inclusive-scans by adding theresult of the left neighbors into the local inclusive-scans. The
last addition can be skipped if m ≤ w − α + 1.
this idea, the grid and thread block levels are omitted in thediscussion below. We consider here cascading the disjoint
segmented scheme with the sequential bit-parallel scheme,
each running at the warp and thread levels, respectively. The
cascaded algorithm executes in the following three phases.
1)

2)

3)

Local inclusive-scan per thread. The l-th thread of
the i-th warp, where 1 ≤ l ≤ 32 and 1 ≤ i ≤⌈n/sα⌉, sequentially computes A′s(i−1)+l , which isthe sum of α contiguous elements of a segment,
def ∑αwhere A′j =l=1 Aα(j−1)+l . Note that each warp
is responsible for sα elements (instead of s elements) owing to algorithm cascading.
Parallel inclusive-scan per warp. According to Algorithm 3, the i-th warp, where 1 ≤ i ≤⌈n/sα⌉ performs the parallel inclusive-scan algorithm of A′s(i−1)+1 , A′s(i−1)+2 , · · · , A′s(i−1)+s to generate every α elements within a segment, i.e.,
′
′
′
Ci,1
, Ci,2, · · · , Ci,s. Note that these s elements arepart of the entire segment; the remaining s(α − 1)elements still hold local sums.
Local addition per thread. To obtain the remaining
global sums, the l-th thread of the i-th warp, where1 ≤ l ≤ 32 and 1 ≤ i ≤ ⌈n/sα⌉, updates its
local sums by adding the global sum of its left
neighboring thread; the global sum to be added is′given by Ci,l−1, if l ≥ 2, otherwise it is given by′
Ci−1,s .

The last addition phase can be skipped if m ≤ w − α + 1.
The key idea behind this elimination is to extend the patternsuch that it ends with α − 1 wildcards, which match toarbitrary characters of the text. Obviously, this extension
reduces the maximum length of the pattern from w to
w − α + 1; however, as shown in Fig. 5, wildcards allow
bit-vector Rj of Eq. (1) to save the current automaton states
for the future α − 1 steps by extending the bit-vector with…

…

…
R j ✏ ( s ✎1)✍R j ✌☞

Rj✡

☛…

m ☎1 m

✝ ✆1

…

✦

j ✄1

…

Rj
m ✟ ✠ ✞1+

Ci4, s✘Ci ,( s ★1)✩ ✧Thread1…

✡

Thread 4 Ci✤, s

+

…
…R

…

+

…

R j ✁✂ 1+…

Ci✣✢1, 2++…

+As✒(i ✑ 2 )✏ 2

…

…

+

+

+

s

…

s ( i 1) s

+

…

+

+

Phase 1 and 2 of Fig. 3… A☞ ☛ ✡ … A✞ ✠ ✟ ✞s ( i 1) 1

9

2 1

m

Fig. 5: Wildcard extension for skipping the last addition
phase. After the partial inclusive-scan algorithm of Fig. 4
executes, Rj holds global results, whereas its preceding
α − 1 elements Rj+α−1 , Rj+α−2 , · · · , Rj−1 hold local results. The global results of the α − 1 elements can be
estimated from those of Rj by extending the vector by α − 1
bits.

α − 1 bits, i.e., for all 0 ≤ i < α and α < j ≤ n,
Rj [m + i] = Rj−i [m]. This equivalent relation is usefulfor retrieving the entire segment from every α elementthat hold the complete inclusive-scans. The j -th thread isallowed to substitute Rj [m], Rj [m + 1], . . . , Rj [m + α − 1]
for Rj [m], Rj+1 [m], . . . , Rα−1 [m], respectively.
As for the appropriate value of α, we experimentally determined the best value that maximized search throughput.
As a guideline, α is at least four for 1-byte data. One may
feel that α > 4 fails to coalesce memory transactions because
a contiguous memory region of larger than 128 bytes issimultaneously accessed by a warp; however, our preliminary results indicate that the number of global memorytransactions does not increase if the memory region size is
512 bytes, which may be due to the demand from computer
graphics applications, which usually deploy the float4
data structure to pack RGBA data into 512-byte data.
5.2

Complete Loop Unrolling

Our proposed method has two distinct loops within the kernel function. One is the d-loop at line 4 of Algorithm 2; theother is the j -loop at line 8, required for the inclusive-scanoperation. As mentioned in Section 4.3, every thread iterates
through the latter loop log m times instead of log n times,
while the former loop executes k times. Search throughput
can be increased by completely unrolling these loops. To
do so, we use the C++ template presented in [40]. Morespecifically, the unrolled kernel function for specific valuesof m and k is automatically generated by the C++ template at compile time. This template-based scheme generatesm × (k + 1) kernel functions. Because constant values of mand k are embedded at compile time, the generated kernel
consumes fewer register files.
One drawback of the template-based scheme here isthat it is not robust against increases in m and k . Let M
and K be the maximum values of m and k acceptable
at compile time, respectively, with m ≤ M and k ≤ K .
Consider a kernel function generated for pattern size m andmaximum error k . The generated function then contains

=32
=16
=8
=4
=2

0

8

16

24
32
40
Pattern size m

48

56

64

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

10
=32
=16
=8
=4
=2

0

8

16

24
32
40
Pattern size m

(a)

48

56

Throughput (Tbps)

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

Throughput (Tbps)

Throughput (Tbps)

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

64

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

=32
=16
=8
=4
=2

0

8

(b)

16

24
32
40
Pattern size m

48

56

64

(c)

Fig. 6: Search throughput of our proposed cuShiftOr implementation with a variety of values for granularity size α. Results
are shown for (a) k = 0, (b) k = 1, and (c) k = 2. Text and alphabet sizes were fixed to n = 231 and σ = 64, respectively.

k log m statements because there is a double-nested loop,
where the outer and inner loops iterate k and log m times,
respectively. Hence, the total number of statements is given∑K ∑log M
by k=1 m=1 k log m. Thus, the code size is bounded byO(K 2 M log M ), which restricts m and k .

6

E XPERIMENTAL R ESULTS

We evaluated our proposed cuShiftOr implementation interms of search throughput ρ. Note that search throughput
does not include data transfer time between the CPU andGPU. Our experimental machine had an Intel Xeon CPU
E5-2660 v3 CPU with 64 GB RAM and an NVIDIA GeForceGTX TITAN X (Maxwell) GPU with 12 GB VRAM. We used
CUDA 7.5 [20], GPU Driver 353.82, and Visual Studio 2013
running on Windows 7.
For a long string that exceeds the capacity of GPU memory, ρ is obviously limited by the bandwidth of the PCIe bus.
On the machine used in our experiments, the data transfer
rate from the CPU to GPU was 10.6 GB/s and 11.2 GB/s forthe opposite direction. By contrast, the effective bandwidthof the GPU memory reached 248.8 GB/s on the experimental
GPU according to the bandwidthTest program [39]. Because
string matching is usually a memory-intensive operation,
we focused on the in-core situation rather than the outof-core situation. Therefore, in-core search throughput ρ isbounded by 248.8 GB/s or 1.99 Tbps.
As for a comparable method, we implemented GBPR
[11] as a segmentation-based method. Our method differsfrom GBPR in realizing (1) a parallel scan algorithm thatremoves duplicate searches by eliminating the halo regionsand (2) algorithm cascading that facilitates coalesced memory access by mapping appropriate algorithms to the thread
hierarchy. We implemented GBPR by ourselves because it
was not publicly available. Our GBPR implementation used
the same data structure as the original GBPR; we usedshared memory to store bit masks and global memory tostore the text and search results. Segment size s was experimentally determined to be s = 8, which maximized searchthroughput ρ on our experimental machine. The readers canrefer to the supplementary material for detailed description
on GBPR including a sensitivity study on s and efficiency
analysis. We also implemented an extended version of GBPR
that realizes coalesced memory access with a tiling technique. Similar to cuShiftOr, the tiled GBPR allows every
thread to fetch α (> 1) contiguous characters at a time tolocally update automata states α times. This procedure isrepeated until the end of its responsible segment is reached.

6.1

Parameter Tuning

We first conducted a preliminary experiment to determine
the best task granularity α that maximized the searchthroughput ρ for cuShiftOr. Figure 6 shows search throughput ρ measured while varying task granularity α andmaximum error k , with maximum text size n = 231 andmaximum alphabet size σ = 64. The text and patterns
were generated randomly. As shown in the figure, α greatly
impacts search throughput ρ, which ranged from 0.26 to
1.88 Tbps when k = 0 (i.e., exact matching). According tothese results, we decided to use α = 16, which achieved thehighest throughput for most pattern sizes.
In Fig. 6a, we observe that cuShiftOr significantly
dropped search throughput to 1.05 Tbps when m > 32. For
such a long pattern, cuShiftOr uses a 64-bit word to store the
bit-vector but the current CUDA does not provide shuffle
instructions for 64-bit data. Consequently, switching word
size w degraded the efficiency of the shuffle instructions
required for inclusive-scan computation; in other words,
a 64-bit shuffle operation was implemented using 32-bit
shuffle instructions. This performance degradation was notobserved with smaller α, where search throughput ρ wasdetermined by off-chip memory access rather than instruction issue rate (i.e., on-chip memory access). Therefore, we
conclude here that choosing the best task granularity α is
key to maximizing search throughput ρ.
6.2

Performance Analysis

Using task granularity α determined using the methods
above, we measured search throughput ρ while varying text
size n, pattern size m, maximum error k , and alphabet size
σ . The text and patterns were again generated randomly.
Figure 7 shows search throughput ρ measured whilevarying pattern size m and maximum error k with maximum text size n = 231 and maximum alphabet size σ = 64.
Here, cuShiftOr was 16.9 times faster than GBPR (α = 1)
when k = 0. In particular, search throughput ρ for a shortpattern (i.e., k = 0 and m ≤ 32) reached 1.88 Tbps, demonstrating a high efficiency of 94% in terms of off-chip memory
R/W throughput. This achieved efficiency indicates that
cuShiftOr successfully exploits the strength of the latency
hiding architecture; in other words, scan computation is
fully overlapped with off-chip memory access.
Both cuShiftOr and GBPR saw degradation in searchthroughput ρ as we increased pattern size m from 33 to 64;
however, this decrease was relatively small for cuShiftOr.
As an example, when k = 0, cuShiftOr decreased search

cuShiftOr
GBPR ( =8)
GBPR ( =4)
GBPR ( =2)
GBPR ( =1)

0

8

16

24 32 40
Pattern size m

48

56

64

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

11

cuShiftOr
GBPR ( =8)
GBPR ( =4)
GBPR ( =2)
GBPR ( =1)

0

8

16

24 32 40
Pattern size m

(a)

48

56

Throughput (Tbps)

2.2
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

Throughput (Tbps)

Throughput (Tbps)

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

64

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

cuShiftOr
GBPR ( =8)
GBPR ( =4)
GBPR ( =2)
GBPR ( =1)

0

8

(b)

16

24 32 40
Pattern size m

48

56

64

(c)

0.01
0.001
cuShiftOr
GBPR ( =1)0.00010.00001

1
0.1
0.01
0.001
cuShiftOr
GBPR ( =1)0.0001

1
0.1
0.01
0.001
cuShiftOr
GBPR ( =1)0.00010.00001

0.00001
210 212 214 216 218 220 222 224 226 228 230 232
Text size nThroughput (Tbps)

1
0.1

Throughput (Tbps)

Throughput (Tbps)

Fig. 7: Search throughput with different pattern sizes m. Results are shown for (a) k = 0, (b) k = 1, and (c) k = 2. Text andalphabet sizes were fixed to n = 231 and σ = 64, respectively. The original GBPR corresponds to α = 1.

210 212 214 216 218 220 222 224 226 228 230 232
Text size n

(a)

(b)

210 212 214 216 218 220 222 224 226 228 230 232
Text size n(c)Fig. 8: Search throughput with varying text sizes n. Results are shown for (a) k = 0, (b) k = 1, and (c) k = 2. Pattern sizem and alphabet size σ were fixed to m = 64 and σ = 64, respectively. The vertical axis is presented as a logarithmic scale.
throughput ρ by only 2%, whereas GBPR (α = 1) decreased ρ by 45%. We attribute this substantial difference
to the deployed parallel scheme. As mentioned above, theinstruction issue rate limited the performance of cuShiftOr.
The number of instructions per thread is given by O(log m)for cuShiftOr because threads compute inclusive-scans inlog m parallel steps. By contrast, threads in GBPR (i.e., thesegmentation-based scheme) are responsible for s + m − 1characters in the text, meaning that each thread must reada segment of O(s + m) characters from off-chip memory,where s = 8. Therefore, cuShiftOr is more robust than
GBPR against increases of pattern size m. Obviously, thisconclusion assumes that pattern size is short enough to
fit within machine word size w. Note that cuShiftOr also
deploys the segmentation-based scheme, but segments are
assigned to thread blocks, namely, a higher level of thread
hierarchy. Therefore, s is much larger than m, which can beignored at the thread level.
As for approximate string matching, cuShiftOr dropped
search throughput ρ as we increased maximum error k .
For example, the maximum ρ for k = 1 (or k = 2) wasapproximately 39% (or 29%) lower than that of k = 0 (ork = 1, respectively). By contrast, GBPR (α = 1) maintained
a high search throughput ρ for arbitrary k . We attribute
this difference to the performance bottleneck of the GPU
code. The performance bottleneck of cuShiftOr was onchip memory, whereas that of GBPR was off-chip memory.
The amount of on-chip memory access increases with kbecause the time complexity of the WM algorithm is givenby O(nm/wk), which linearly increases with k . By contrast,
the amount of off-chip memory access is independent ofk because Algorithm 2 reads and writes off-chip memory
data only once at lines 1 and 13, respectively. Consequently,
cuShiftOr degraded search throughput ρ as k increased,
whereas GBPR held ρ almost constant for arbitrary k .
Our tiled version of GBPR (α > 1) explains why

cuShiftOr was faster than the original GBPR. In Fig. 7c, thetiled GBPR (α = 8) achieved up to 93% of search throughput
obtained via cuShiftOr; as for approximate string matching,
the tiling technique (i.e., memory coalescing) is key for eliminating the performance gap between cuShiftOr and GBPR.
However, memory coalescing was insufficient in yielding
Tbps speeds for search throughput where on-chip memory
access determines performance. In this case, the amountof on-chip memory access must be reduced by eliminating
duplicate searches.
Figure 8 shows search throughput measured while varying text size n for maximum pattern size m = 64 andσ = 64. The figure indicates that the text should contain atleast 224 characters (i.e., 16 MB text data) to achieve at least
75% of the maximum search throughput. In contrast, GBPR
achieved 75% of the maximum search throughput with 222
characters, which was 25% smaller than that needed for
cuShiftOr. This difference occurred due to a side effect ofalgorithm cascading in which a large α causes serialization
on each thread. In fact, the ratio for 222 characters wasincreased from 49% to 74% by reducing task granularity α
from 16 to 2. Therefore, cuShiftOr requires a larger text than
GBPR to increase the ratio of parallel processing over serial
processing, as Amdahl’s law [41] points out.
We also found that, when n ≤ 222 , GBPR varied ρ
according to n, which differed from the behavior that wetheoretically expected. GBPR processes ⌈n/s⌉ segments ofsize s + m − 1 on p processors. Consequently, assumingthat m = O(w) (i.e., short patterns), the search throughput of GBPR is expected to be sp/(s + m − 1), which is
independent of n. This gap between theoretical behavior
and actual behavior is due to the lack of parallelism thatcan cause idle time on massively parallel GPU threads.
The GPU architecture requires a large number of threads
to sufficiently hide off-chip memory latency by overlapping
memory access with data-independent computation. NoticecuShiftOr
GBPR (✁=1)

0

64

128
Alphabet size

192

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

256

12

cuShiftOr
GBPR (✁=1)

0

64

128
Alphabet size

(a)

192

Throughput (Tbps)

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

Throughput (Tbps)

Throughput (Tbps)

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

256

cuShiftOr
GBPR (✁=1)

0

64

(b)

128
Alphabet size

192

256

(c)

6.3

Performance Comparison with Corpus

In this subsection, we show performance comparison results
obtained using practical data sets. For text data, we used the
Pizza&Chilli Corpus [42], which includes corpus data sets
across a variety of fields. The first 64 characters of the text
were selected as the pattern to be searched for. In addition
to GBPR [11], we also compared cuShiftOr with XBitPar [26]obtained from http://xbitpar.sourceforge.net/.
Figure 10 shows measured search throughputs. Whenk = 0 and m = 64, search throughput ρ of cuShiftOr
ranged from 1782 to 1872 Gbps. These throughputs were
30% higher than those expected from results obtained withrandomly-generated data sets (i.e., Fig. 9a); we expected
search throughput to be approximately 1400 Gbps whenσ > 32. This unexpected behavior was due to the skewed
distribution of characters in the corpus, which did notexist in the randomly-generated data sets. As an example,
the ENGLISH data set contained 239 characters, but its
inverse probability of matching was 28.7, indicating thatspecific characters in the alphabet frequently appeared in
the text. In such a case, cuShiftOr can slightly increase searchthroughput by avoiding bank conflicts in shared memory.
For example, bit masks for A, G, T, and C are stored in
the 1st, 7th, 20th, and 3rd banks, respectively, because thetable of bit masks is indexed using the corresponding ASCII
codes. Consequently, bank conflicts will not occur for biological data that contains these four characters. In contrast,
randomly generated datasets contain 256 characters, thus
each bank is responsible for eight characters. In this case,bank conflicts occur when 32 threads in a warp access thesecharacters at the same time.

10000

cuShiftOr

GBPR

XBitPar

1000
100
10
1
4

8 16 32 64 4
PITCHES8 16 32 64 4 8 16 32 64 4 8 16 32 64 4 8 16 32 64 4
XML
SOURCESPROTEINSENGLISH
Pattern size m

8 16 32 64
DNA

Throughput (Gbps)

(a)
10000

cuShiftOr

GBPR

XBitPar

1000
100
10
1
4

8 16 32 64 4
PITCHES8 16 32 64 4 8 16 32 64 4 8 16 32 64 4 8 16 32 64 4
XML
SOURCESPROTEINSENGLISH
Pattern size m

8 16 32 64
DNA(b)Throughput (Gbps)

that cuShiftOr also faced on the same issue when n ≤ 222 .
Finally, as shown in Fig. 9, we measured search throughput ρ while varying alphabet size σ for maximum text sizen = 231 and maximum pattern size m = 64. When k = 0,
ρ clearly decreased as we increased σ . By contrast, we did
not observe such behavior when k > 0. Search throughput
ρ for k = 0 was determined by on-chip memory accessrather than off-chip memory access. Consequently, bankconflicts in shared memory limited ρ when σ > 32 because
the shared memory of the current CUDA is structured
with 32 banks; further, the cuShiftOr implementation uses
shared memory to store a table indexed by alphabet Σ,as mentioned in Section 5. Realizing conflict-free access tothis shared table is not easy owing to the irregularity ofcharacters in the text.

Throughput (Gbps)

Fig. 9: Search throughput with varying alphabet sizes σ . Results are shown for (a) k = 0, (b) k = 1, and (c) k = 2. Text sizeand pattern size were fixed to n = 231 and m = 64, respectively.

10000

cuShiftOr

GBPR

XBitPar

1000
100
10
1
4

8 16 32 64 4
PITCHES8 16 32 64 4

XML

8 16 32 64 4 8 16 32 64 4 8 16 32 64 4
SOURCESPROTEINSENGLISH
Pattern size m

8 16 32 64

DNA

(c)

Fig. 10: Search throughput using the Pizza&Chilli Corpus[42]. Results are shown for (a) k = 0, (b) k = 1, and (c)k = 2, where 4 ≤ m ≤ 64. PITCHES, XML, SOURCES,PROTEINS, ENGLISH, and DNA data sets have alphabet
sizes σ of 133, 97, 230, 27, 239, and 16, respectively, and theirinverse probabilities of matching are 39.8, 28.7, 24.8, 17.0,
15.3, and, 3.9, respectively.

The increases in terms of runtimes of cuShiftOr over
GBPR were similar to those obtained with the random data
presented in Section 6.2 above. Here, cuShiftOr was 15.7–
16.7 times faster than GBPR when m = 32 and k = 0, but
the gap between performance measures decreased as we
increased k . More specifically, we observed 9.3–10.2 times
and 6.7–7.4 times when k = 1 and k = 2, respectively. Thus,
the family of bit-parallel algorithms demonstrated stablesearch throughput for arbitrary data.
Both cuShiftOr and GBPR failed to search a pattern
longer than 64 characters, whereas XBitPar successfully
processed a long pattern of up to 1024 characters. Nonetheless, cuShiftOr was 93–108 times faster than XBitPar whenm = 32 and k = 0. XBitPar adopts a 1024-bit virtual
word for arbitrary m, thus the execution efficiency drops for
short queries. More specifically, w − m bits in a word were
wasted during search, thus only 3% of threads were utilized
when m = 32. Therefore, cuShiftOr and XBitPar should beappropriately switched according to pattern length m. Note

6.4

Performance Comparison with Sequence Aligners

Finally, we compared cuShiftOr with indexing-based sequence aligners running on a CPU: BitMapper [34] andmrFAST [35] as all-mappers, which find all mappings of agiven sequencing read in the reference genome by setting
a maximum edit distance, and CUSHAW2 [36], BWA-MEM[37], and Bowtie2 [38] as best-mappers, which aim to report
the best few mappings under similar constraints. Millionsof short reads from the 1000 Genomes Project were aligned
to Homo sapiens chromosome 1, GRCh38.p7. As shownin Fig. 11, cuShiftOr achieved search throughputs of 0.35–
0.4 Tbps. In contrast, sophisticated indexes allowed thealigners to focus on possible matching locations in the text.
Consequently, all-mappers and best-mappers aligned reads
at approximately 3–12 Tbps and 132–836 Tbps, respectively.
Figure 11 also shows mapping ratios, i.e., the numberof mapped reads divided by the total number of reads.
The mapping ratios of cuShiftOr were not exactly the sameas those of mrFAST [35], although the same maximum
edit distance of two was commonly used for search. Thisdifference was due to reverse complementary sequences andambiguous nucleotides; cuShiftOr accepts a single pattern
per search, thus the search process must be invoked twice
for the original pattern and its reverse complementary pattern. Moreover, these results must be merged into a single
result. Finally, cuShiftOr neither skips reads that include
many ambiguity characters nor finds a match between the
ambiguity character (N) and any character (A, G, T, or C).
Consequently, we think that cuShiftOr is useful for accelerating online searching where indexes cannot be constructed before searching. For example, cuShiftOr can beused for text stream mining in social networks, where text
streams are analyzed on the fly. Because CPU-GPU datatransfer usually limits the performance of GPU systems, theimpact of cuShiftOr increases if text streams are iteratively
analyzed on the GPU. Similar promising examples are network intrusion detection and code debugging, where log
data are iteratively examined to find regions of interest.

7

C ONCLUSION

In this paper, we presented a tribrid parallel method, named
cuShiftOr, that is capable of exploiting massive data parallelism inherent in string matching. Our proposed methodinterprets the SO and WM algorithms as inclusive-scanoperations. This interpretation allows massive numbers of

258.930.40

N/A

284.58

3.39

20%
18%
16%
14%
12%
10%
8%
6%
4%
2%
0%

1000
900
800
700
600
500
400
300
200
100
0

639.66218.69132.23
0.35

10.70

12.39

cuShiftOr BitMapper mrFast CUSHAW2 BWA Bowtie 2

cuShiftOr BitMapper mrFast CUSHAW2 BWA Bowtie 2

(a)

(b)

20%
18%
16%
14%
12%
10%
8%
6%
4%
2%
0%

Mapping ratio

835.85

Search throughput (Tbps)

1000
900
800
700
600
500
400
300
200
100
0

Search throughput (Tbps)

that cuShiftOr can be extended to search a long pattern witha virtual word. However, this extension degrades searchthroughput due to the overhead required to move data
between physical words.
Thus, the limitation on pattern length is critical for long
patterns, but there are many real-world cases that require
support only for short patterns. For example, a text search
query typically consists of two to four words and the patternlength distribution follows a power law [43]. With respect
to biological applications, next-generation sequencing systems produce short reads up to several hundred characters,
which are then aligned to a long sequence of millions ofcharacters.

13

Mapping ratio

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

Fig. 11: Alignment results of cuShiftOr and sequence aligners [34]–[38]. Results for (a) National Center for Biotechnology Information (NCBI) SRR003084 (8.8 M reads of 36 basepair (bp)) and (b) NCBI SRR003092 (16.1 M reads of 51 bp).
BitMapper failed to align SRR003084.
threads to run without halo regions, but requires global synchronization. To reduce these synchronization costs, we integrated the scan scheme into a segmentation-based scheme,
which does require halo regions (i.e., duplicate work) but
avoids synchronization. We showed that these contrasting
methods can be appropriately mapped onto the hierarchyof GPU threads. As such, the scan scheme exploits finegrained parallelism to minimize duplicate work, whereas
the segmentation-based scheme exploits coarse-grained parallelism to avoid global synchronization.
We analyzed the performance of cuShiftOr and compared its performance to previous segmentation-based
methods. In our results, cuShiftOr achieved 1.88 Tbps ofin-core search throughput for short queries, which was
equivalent to 94% execution efficiency in terms of off-chip
memory R/W throughput. However, the throughputs were
not satisfactory compared to those of sequence aligners that
construct indexes before searching. Thus, we conclude that
cuShiftOr is useful for accelerating online string matchingfor short patterns containing up to 64 characters.

ACKNOWLEDGMENTS
This study was supported in part by the Japan Society
for the Promotion of Science KAKENHI Grant Numbers
15K12008, 15H01687, and 16H02801. We are also gratefulto the anonymous reviewers for their valuable comments.

R EFERENCES
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]

A. Tumeo and O. Villa, “Accelerating DNA analysis applications
on GPU clusters,” in Proc. 8th Symp. Appl. Specific Processors, Jun.
2010, pp. 71–76.
S. Faro and T. Lecroq, “The exact online string matching problem:
A review of the most recent results,” ACM Comput. Surv., vol. 45,
no. 2, Feb. 2013, article 13, 42 pages.
C.-H. Lin, C.-H. Liu, L.-S. Chien, and S.-C. Chang, “Accelerating
pattern matching using a novel parallel algorithm on GPUs,” IEEE
Trans. Comput., vol. 62, no. 10, pp. 1906–1916, Oct. 2013,
G. Navarro, “A guided tour to approximate string matching,”
ACM Comput. Surveys, vol. 33, no. 1, pp. 31–88, Mar. 2001.
D. E. Knuth, J. H. Morris, and V. R. Pratt, “Fast pattern matchingin strings,” SIAM J. Comput., vol. 6, no. 2, pp. 323–350, Jun. 1977.
A. V. Aho and M. J. Corasick, “Efficient string matching: An aid tobibliographic search,” Commun. ACM, vol. 18, no. 6, pp. 333–340,
Jun. 1975.
R. Baeza-Yates, “Efficient text searching,” Ph.D. dissertation, University of Waterloo, Waterloo, ON, Canada, May 1989.
R. Baeza-Yates and G. H. Gonnet, “A new approach to text
searching,” Commun. ACM, vol. 35, no. 10, pp. 74–82, Oct. 1992.
S. Wu and U. Manber, “Fast text searching allowing errors,”
Commun. ACM, vol. 35, no. 10, pp. 83–91, Oct. 1992.

IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. XX, XX 201X

[10] K. Kusudo, F. Ino, and K. Hagihara, “A bit-parallel algorithm
for searching multiple patterns with various lengths,” J. Parallel
Distrib. Comput., vol. 76, pp. 49–57, Feb. 2015.
[11] C. H. Lin, G. H. Wang, and C. C. Huang, “Hierarchical parallelism of bit-parallel algorithm for approximate string matching
on GPUs,” in Proc. Symp. Comput. Appl. and Commun., Jul. 2014,
pp. 76–81.
[12] T. T. Tran, S. Schindel, Y. Liu, and B. Schmidt, “Bit-parallel approximate pattern matching on the Xeon Phi coprocessor,” in Proc. 26th
Int. Symp. Comput. Archit. and High Perform. Comput., Oct. 2014, pp.
81–88.
[13] G. E. Blelloch, “Scans as primitive parallel operations,” IEEE Trans.
Comput., vol. 38, no. 11, pp. 1526–1538, Nov. 1989.
[14] M. Harris, S. Sengupta, and J. D. Owens, “Parallel prefix sum(scan) with CUDA,” in GPU Gems 3, H. Nguyen, Ed. Addison
Wesley, Aug. 2007, ch. 39, pp. 851–876.
[15] A. Khajeh-Saeed, S. Poole, and B. Perot, “Acceleration of theSmith-Waterman algorithm using single and multiple graphics
processors,” J. Comput. Physics, vol. 229, no. 11, pp. 4247–4258,
Jun. 2010.
[16] T. F. Smith and M. S. Waterman, “Identification of common molecular subsequences,” J. Molecular Biol., vol. 147, no. 1, pp. 195–197,
Mar. 1981.
[17] S. Sengupta, M. Harris, Y. Zhang, and J. D. Owens, “Scan primitives for GPU computing,” in Proc. 22nd ACM SIGGRAPH/EUROGRAPHICS Symp. Graph. Hardware, Aug. 2007, pp. 97–106.
[18] G. E. Blelloch, “Prefix sums and their applications,” School ofComputer Science, Carnegie Mellon University, Tech. Rep. CMUCS-90-190, Nov. 1990.
[19] W. D. Hillis and G. L. Steele, “Data parallel algorithms,” Commun.
ACM, vol. 29, no. 12, pp. 1170–1183, Dec. 1986.
[20] NVIDIA Corporation, “CUDA C Programming Guide Version
7.5,” Sep. 2015. [Online]. Available: http://docs.nvidia.com/
cuda/pdf/CUDA C Programming Guide.pdf[21] ——, “NVIDIA GeForce GTX 980,” Nov. 2014. [Online]. Available:
http://international.download.nvidia.com/geforce-com/
international/pdfs/GeForce GTX 980 Whitepaper FINAL.PDF[22] Y. Liu, L. Guo, J. Li, M. Ren, and K. Li, “Parallel algorithms forapproximate string matching with k mismatches on CUDA,” inProc. IEEE 26th Int. Symp. Parallel and Distrib. Process. Workshops
and PhD Forum, May 2012, pp. 2414–2422.
[23] D. Man, K. Nakano, and Y. Ito, “The approximate string matchingon the hierarchical memory machine, with performance evaluation,” in Proc. IEEE 7th Int. Symp. Embedded Multicore/Manycore
System-on-Chip, Sep. 2013, pp. 79–84.
[24] K. Xu, W. Cui, Y. Hu, and L. Guo, “Bit-parallel multiple approximate string matching based on GPU,” in Proc. 1st Int. Conf. Inf.
Technol. and Quantitative Manage., May 2013, pp. 523–529.
[25] M. Onsjö and O. Watanabe, “Implementation of a bit-parallel
approximate string matching algorithm,” in IPSJ SIG Technical
Report, vol. 2009-AL-124, no. 3, May 2009.
[26] T. T. Tran, Y. Liu, and B. Schmidt, “Bit-parallel approximate patternmatching: Kepler GPU versus Xeon Phi,” Parallel Comput., vol. 54,
pp. 128–138, May 2016.
[27] Intel
Corporation,“Intelarchitecture
instruction
setextensions programming reference,” Dec. 2013. [Online].
Available:
http://download-software.intel.com/sites/default/
files/managed/71/2e/319433-017.pdf[28] R. Prasad, S. Agarwal, I. Yadav, and B. Singh, “A fast bitparallel multi-patterns string matching algorithm for biologicalsequences,” in Proc. Int. Symp. Biocomput., no. 46, Feb. 2010, 4
pages.
[29] R. Baeza-Yates and G. Navarro, “Faster approximate string matching,” Algorithmica, vol. 23, no. 2, pp. 127–158, Feb. 1999.
[30] G. Myers, “A fast bit-vector algorithm for approximate string
matching based on dynamic programming,” J. ACM, vol. 46, no. 3,
pp. 395–415, May 1999.
[31] L. Langner, “Parallelization of Myers fast bit-vector algorithm
using GPGPU,” Ph.D. dissertation, Freie Universität Berlin, Berlin,
Germany, Apr. 2011.
[32] D. Man, K. Nakano, and Y. Ito, “An optimal implementationof the approximate string matching on the hierarchical memorymachine, with performance evaluation on the GPU,” IEICE Trans.
Inf. Syst., vol. E97-D, no. 12, pp. 3063–3071, Dec. 2014.
[33] C.-H. Lin, S.-Y. Tsai, C.-H. Liu, S.-C. Chang, and J.-M. Shyu,“Accelerating string matching using multi-threaded algorithm on
GPU,” in Proc. Global Commun. Conf., Dec. 2010, 5 pages.

14[34] H. Cheng, H. Jiang, J. Yang, Y. Xu, and Y. Shang, “Bitmapper:
an efficient all-mapper based on bit-vector computing,” BMC
Bioinformatics, vol. 16, no. 192, Jun. 2015.
[35] C. Alkan, J. M. Kidd, T. Marques-Bonet, G. Aksay, F. Antonacci,
F. Hormozdiari, J. O. Kitzman, C. Baker, M. Malig, O. Mutlu,
S. C. Sahinalp, R. A. Gibbs, and E. E. Eichler, “Personalized copy
number and segmental duplication maps using next-generation
sequencing,” Nature Genetics, vol. 41, no. 10, pp. 1061–1067, Aug.
2009.
[36] Y. Liu and B. Schmidt, “Long read alignment based on maximum
exact match seeds,” Bioinformatics, vol. 28, no. 18, pp. i318–i324,
Sep. 2012.
[37] H. Li, “Aligning sequence reads, clone sequences and assembly
contigs with BWA-MEM,” arXiv:1303.3997 [q-bio.GN], 2013.
[38] B. Langmead and S. L. Salzberg, “Fast gapped-read alignment
with Bowtie 2,” Nature Methods, vol. 9, no. 4, pp. 357–359, Mar.
2012.
[39] NVIDIA Corporation, “CUDA Code Samples,” 2014. [Online].
Available: http://developer.nvidia.com/cuda-code-samples/.
[40] M. Harris, “Optimizing parallel reduction in CUDA,” Nov.
2007. [Online]. Available: http://developer.download.nvidia.
com/assets/cuda/files/reduction.pdf.
[41] G. M. Amdahl, “Validity of the single processor approach toachieving large scale computing capabilities,” in Proc. the AFIPS
Conf., 1967, pp. 483–485.
[42] P. Ferragina and G. Navarro, “Pizza&chili corpus compressedindexes and their testbeds,” Sep. 2005. [Online]. Available: http://pizzachili.dcc.uchile.cl/.
[43] A. Arampatzis and J. Kamps, “A study of query length,” inProc. 31st Ann. Int. ACM SIGIR Conf. Res. and Development in Inf.
Retrieval, Jul. 2008, pp. 811–812.

PLACE
PHOTO
HERE

Yasuaki Mitani received the B.E. and M.E.
degrees in information and computer sciencesfrom Osaka University, Osaka, Japan, in 2014
and 2016, respectively. He is currently with
DWANGO Corporation, Ltd. His current researchinterests include high performance computing
and string processing.

PLACE
PHOTO
HERE

Fumihiko Ino (S’01–A’03–M’04) received the
B.E., M.E., and Ph.D. degrees in informationand computer sciences from Osaka University,
Osaka, Japan, in 1998, 2000, and 2004, respectively. He is currently an Associate Professor in
the Graduate School of Information Science andTechnology at Osaka University. His research interests include parallel and distributed systems,
software development tools, and performance
evaluation.

PLACE
PHOTO
HERE

Kenichi Hagihara received the B.E., M.E., andPh.D. degrees in information and computer sciences from Osaka University, Osaka, Japan, in
1974, 1976, and 1979, respectively. From 1994
to 2002, he was a Professor in the Graduate
School of Engineering Science, Osaka University. Since 2002, he has been a Professor in
the Graduate School of Information Science andTechnology, Osaka University. His research interests include the fundamentals and practicalapplication of parallel processing.