Accelerating Scoring Computation of Smith-Waterman Algorithm with
Mixed Word Length
Kazuki Yasui and Fumihiko Ino
Graduate School of Information Science and Technology, Osaka University
1-5 Yamadaoka, Suita, Osaka 565-0871, Japan
Email: {k-yasui, ino}@ist.osaka-u.ac.jp
Abstract—In this paper, we propose a graphics processing unit
(GPU) accelerated Smith-Waterman (SW) method for aligninga short sequence with a long sequence. Our method deploys
an 8-bit data structure to accelerate scoring computation,
which can be classified into a memory-intensive operation.
The proposed method is based on a mixed-word-length scheme
capable of appropriately switching 8-bit mode and 32-bit mode
to avoid integer overflow during computation. We integrate
our data structure and scheme into a scan-based parallel
approach to achieve high throughput for short sequences,which usually limit the parallelism inherent in the SW computation. Experimental results show that the proposed method
is approximately three times faster than a previous scanbased approach that uniformly deploys 32-bit data structure.
We also show that, for short sequences, scan-based parallel
approaches are about twenty times faster than an antidiagonalbased parallel approach typically deployed for pairs of longsequences.

1. Introduction
Pairwise sequence alignment is a process that identifies
similar regions between two biological sequences such as
nucleotide and protein sequences. This fundamental process is useful when analyzing functional, structural, andevolutional relationships between two sequences. There are
two types of sequence alignment: global and local alignment. The former produces an end-to-end alignment of
the sequences and the latter produces local regions (i.e.,
subsequences) with the highest similarity. Local alignment
is typically used for constructing an evolutionary tree because this technique identifies regions where mutation like
insertions or deletions of nucleotides occurred in the processof evolution.
The Smith-Waterman (SW) algorithm [1] is commonlyused for finding the optimal local alignment. Given two
sequences of length m and n (≤ m), the SW algorithm
runs in O(mn) time. The most time-consuming part of theSW algorithm is a matrix-filling phase, in which similarity
scores of arbitrary regions of the sequences are computed
as elements of a 2-D scoring matrix. Consequently, acceleration methods are required to apply the SW algorithm togrowing biological sequences, which can have lengths ofgigabase pairs (Gbp).
Owing to this growing demand, some researchers implemented the SW algorithm on accelerators such as thegraphics processing unit (GPU) [2], [3], [4], [5] and Xeon
Phi [6]. Among these promising accelerators, the GPU is acommodity device that has been accelerated memory- andcompute-intensive applications in various fields [7], [8], [9].
A C-like programming framework called compute unified
device architecture (CUDA) [10] is available for application
developers to offload their time-consuming tasks from the
CPU to the NVIDIA GPU.
As for the acceleration of local alignment, SW# [3] and
CUDAlign [4] are useful for massively long sequences of
up to 180 megabase pairs (Mbp) × 183 Mbp. These previous methods are based on an antidiagonal-based parallel
approach [11], [12] in which min(m, n) matrix cells onan antidiagonal are computed in parallel. This approach isuseful for sequences with similar lengths, i.e., m ≃ n. However, the exploited parallelism is significantly limited forshort sequences, for example, 200-bp sequences enormously
produced by short-read sequencing instruments. Achieving
acceleration for such short sequences is a challenging issue
for high-performance computing community because (1) theavailable data parallelism is limited due to small data and(2) recent hardware advances rely on increasing the numberof cores rather than the frequency of clocks. A new strategy
is required to deal with strong scaling cases, which consider
fixed data size on variable numbers of processing elements.
In this paper, we present a GPU-accelerated SW method
for aligning a short sequence with a long sequence. Our
method deploys an 8-bit data structure to save the amountof memory access for matrix filling of the SW algorithm.
Because the 8-bit data structure can suffer from integer overflow, we deploy a mixed-word-length scheme that appropriately selects either 8-bit mode or 32-bit mode according toruntime situation. We also integrate this idea into a scanbased parallel approach [13] in which max(m, n) matrix
cells on a row are computed in parallel. The scan operation
is further localized for optimization.
The rest of the paper is organized as follows. Section 2
gives an overview of the SW algorithm and the scan-based
parallel approach [13]. Section 3 then describes our method
and Section 4 shows experimental results. Finally, Section

5 concludes the paper with future work.

T G T C G

A

T

0

0

0

0

0

0

0

0

C

0

0

0

0

2

1

0

0

T

0

2

1

2

1

1

0

2

G

0

1

4

3

2

3

2

1

T

0

2

3

6

5

4

3

4

A

0

1

2

5

5

4

6

5

C

0

0

1

4

7

6

5

5

2. GPU-Accelerated Pairwise Alignment
The SW algorithm [1] is based on a dynamic programming approach that obtains the optimal local alignment
between two sequences. Figure 1 illustrates how the SW
algorithm computes the local alignment. This algorithm
finds similar subsequences in two phases: (1) computation of a scoring matrix H and (2) backtracking from thematrix element with the highest score. As compared withthe backtracking phase, the matrix-filling phase involves
computation by more than an order of magnitude. Therefore,
we focus on the matrix-filling phase, which we parallelize on
the GPU. The backtracking phase can be rapidly processedon the CPU [14].
Let A = a1 a2 . . . an denote a sequence of length n.
Let B = b1 b2 . . . bm denote a sequence of length m (≥ n)
to be compared with the sequence A. The SW algorithm
then computes an n × m matrix H to obtain the similarity
for any pair of subsequences. Let Ei,j and Fi,j , where1 ≤ i ≤ n and 1 ≤ j ≤ m, be the maximum similarities
involving the first i-th symbols in A and the first j -th
symbols in B , respectively. The maximum similarity Hi,jof two subsequences ending in ai and bj , respectively, isthen recursively defined by
Hi,j = max(0, Ei,j , Fi,j , Hi−1,j−1 + S(ai , bj )),(1)

where S(ai , bj ) represents the score of aligning ai and bj .
The score is typically given as a substitution matrix S . Here,
we define S(ai , bj ) as follows.
{
Smatch ,if ai = bj ,S(ai , bj ) =(2)
Sunmatch , otherwise.

= max(Hi,j−1 − Gs , Ei,j−1 − Ge ),= max(Hi−1,j − Gs , Fi−1,j − Ge ),steps. Therefore, this approach is suited for a pair of shortand long sequences, i.e., m > n.
The scan-based approach fills the matrix in three steps
by rewriting the original recurrences of Eqs. (1)–(4) inthree equations. The first step partially computes the scoring
matrix element H̃i,j as follows.
H̃i,j= max(Hi−1,j−1 + S(ai , bj ), Fi,j , 0),

(3)
(4)

where Gs and Ge are penalties for opening a new gap andfor extending an existing gap, respectively. The values for
Hi,j , Ei,j , and Fi,j are defined as zero if i < 1 or j < 1.
The penalties are called linear gap penalties if Gs = Ge .
Otherwise, they are affine gap penalties. Our method deals
with affine gap penalties. From biological point of view,
staring a gap is more important than the extension of gaps,
i.e., Gs > Ge .

2.1. Scan-based Parallelization
Khajeh-Saeed et al. [13] proposed a scan-based parallel
approach for the SW algorithm. Their approach rewrites
the recurrence relation such that matrix elements on thesame row can be computed in parallel. As compared to theantidiagonal-based approach [11], [12], the exploited parallelism increases from min(m, n) to max(m, n). Similarly,
the number of steps changes from m+n−1 steps to n log m(5)where Fi,j = max(Fi−1,j , Hi−1,j − Gs ) − Ge , 1 ≤ i ≤ n,
and 1 ≤ j ≤ m. Equation (5) indicates that the i-th row
depends only on the (i − 1)-th row. Therefore, this step
can be parallelized by assigning the i-th row to max(m, n)processing elements.
The second step computes the row maximums:
Ẽi,jSimilarities Ei,j and Fi,j are given by
Ei,j
Fi,j

Figure 1. An overview of the SW algorithm. Given two DNA sequences
TGTCGAT and CTGTAC with Smatch = 2 and Sunmatch = −1, thescoring matrix H is computed according to a recurrence relation. After
filling the scoring matrix, the backtracking procedure from the matrix
element with the highest score produces the most similar subsequences:
TGTC and TGTAC.

=

max (H̃i,j−k − k × Ge ),

1<k<j

(6)

where 1 ≤ i ≤ n and 1 ≤ j ≤ m. Equation (6) can beprocessed by performing a scan operation. Given a vectorx = (x1 , x2 , . . . , xn ) of n elements, a scan operation withan associative binary operator ⊕ produces a vector y =(y1 , y2 , . . . , yn ), where yj = x1 ⊕ x2 ⊕ . . . ⊕ xj and 1 ≤j ≤ n. Equation (6) can be represented as a scan operation
with the following associative binary operator ⊕:x ⊕ y = max(x − Ge , y).

(7)

Similar to the first step, the second step can be processedby max(m, n) processing elements in parallel. As for GPUbased implementations, the Thrust library [15] is widelyused to accelerate scan operations on the GPU.
Finally, the scoring matrix Hi,j can be obtained byHi,j=max(H̃i,j , Ẽi,j − Gs ),(8)

where 1 ≤ i ≤ n and 1 ≤ j ≤ m.
The abovementioned three steps must be sequentiallyapplied to every rows from top to bottom of the 2-D scoring
matrix. Only two rows rather than the entire matrix must be
stored in the GPU memory to find the maximum score.

3.2. Mixed-word-length scheme

Hi,0 Hi,1 Hi,2 Hi,3
Pack
0

1

2

…
…

/4 

UnpackHi,0 Hi,1 Hi,2 Hi,3
Figure 2. Proposed packed data structure. m matrix elements are stored
as a 32-bit array of ⌈m/4⌉ elements. Reduced amount of memory accessachieves acceleration although packing and unpacking procedures involve
overheads.

The 8-bit data structure limits the alignment scores in
the range [0, 255]. Consequently, integer overflow occurs if
final scores are greater than 255. Therefore, we have to avoid
integer overflow by appropriately selecting the 8-bit modeor the 32-bit mode according to the runtime situation.
The appropriate mode can be selected for every row
before proceeding to the next row. In more detail, given
the i-th row, the maximum value on the (i + 1)-th row canbe estimated using the substitution matrix S . That is, ourmethod avoids integer overflow by selecting the 32-bit modeif
max Hi,j + Smatch > 255,
(9)
0<j<m

1
2
3
4
5

int x = INT4_PACK(1, 2, 3, 4);
int y = INT4_PACK(5, 6, 7, 8);
int res; // res = x + y = (6, 8, 10, 12)asm("vadd4.u32.u32.u32.sat %0, %1, %2, %3;"
: "=r"(res) : "r"(x), "r"(y), "r"(0));

Figure 3. Example of vector addition of 8-bit data structure. Addition is
implemented with inlining a PTX instruction into the CUDA code.

3. Proposed Method
The proposed alignment method focuses on a pair ofshort and long sequences. We integrate three optimization
strategies into the scan-based parallel approach [13].
1)
2)
3)

8-bit data structure for saving the amount of memory access.
Mixed-word-length scheme for avoiding integer
overflow.
Localized parallel scan for reducing global synchronization.

3.1. 8-bit data structure
Our data structure packs four matrix elements into a
32-bit word on global memory [10] (Fig. 2). Therefore,
the data structure replaces a naive array of m elementswith that of ⌈m/4⌉ elements. This packed data structure
must be referred with 32-bit word access to reduce theamount of memory access by a quarter. In other words, four
consecutive elements on a row are assigned to every thread,allowing threads to access four consecutive elements at atime.
The packed data should be directly processed without any overhead. To realize such direct operations, we
utilize single-instruction, multiple-data (SIMD) instructions
that operate four 8-bit values at a time. Figure 3 shows
an example of vector addition implemented by inlining aparallel thread execution (PTX) assembly statement [10]into the CUDA code. In addition to this addition operator,
the CUDA provides a set of fundamental operators, such as
subtraction and maximization, needed for processing Eqs.
(5)–(8) with the 8-bit data structure.

The overhead for this runtime selection is negligible because
the scan-based parallel scheme computes the maximum
score for every row, as shown in Eq. (8).
On the other hand, our method switches from the 32-bit
mode to the 8-bit mode if
max Hi,j + Smatch ≤ 255 − T,

0<j<m

(10)

where T = 10 is a threshold value determined experimentally. We introduce this threshold to prevent frequentswitches between 8-bit and 32-bit modes because every
switch requires data conversion between 8-bit format and32-bit format. This conversion involves additional overheads, but reducing memory access amount is more important to accelerate memory-intensive operations.
Algorithm 1 shows our mixed-word-length scheme,
which computes the scoring matrix from the top to thebottom rows. As shown in line 4, the scoring matrix is
computed with the appropriate data structure to save theamount of memory access. However, the data conversion atlines 12 and 15 is required to avoid integer overflow. The
highest similarity score and its position are computed atlines 5–7. First, the highest scores within each thread block[10] are computed on the GPU. Secondly, these scores aretransferred to the main memory to reduce them into the finalscore on the CPU.

3.3. Localized parallel scan
The previous method [13] scans the entire data tocompute Eq. (6). However, this global scan can be localized for acceleration because the matrix element Hi,j
is affected by some previous p symbols rather than theentire of max(m, n) = n symbols. Given a vector x =(x1 , x2 , . . . , xn ) of n elements, a local scan operation withan associative binary operator ⊕ produces a vector y =(y1 , y2 , . . . , yn ), where yj = xj−p ⊕ xj−p+1 ⊕ . . . ⊕ xj andp ≤ j ≤ n; the parameter p, the maximum distance that ascore at arbitrary position can affect, is given byp = ⌊n × Smatch /Ge ⌋.

(11)

Equation (11) can be obtained as follows. Equation (7)indicates that the similarity score can decrease by Ge when

Algorithm 1 Mixed-word-length scheme
Input: Sequences A = a1 a2 . . . an , B = b1 b2 . . . bm , and substitution matrix S
Output: Highest similarity score H and its position P
1: mode ← “8BIT”;
2: H ← 0; P ← (0, 0);
3: for k ← 0 to n do
4:Compute the k -th row with mode mode using Eqs. (5)–(8) on the GPU;
5:
Gather the highest scores within each thread block on the GPU;
6:
Transfer the candidate scores from the GPU to the CPU;
7:
Find the highest score Hk and its position Pk ;8:
if Hk > H then9:
H ← Hk ; P ← Pk ;10:
end if11:
if (mode == “8BIT”) and (Eq. (9) is true) then12:
Convert 8-bit data to 32-bit data on the GPU;13:
mode ← “32BIT”;
14:else if (mode == “32BIT”) and (Eq. (10) is true) then15:
Convert 32-bit data to 8-bit data on the GPU;
16:
mode ← “8BIT”;
17:
end if
18: end forprocessing a binary operator ⊕. On the other hand, thehighest possible score on matrix H can reach n × Smatch .
Therefore, this score n × Smatch decreases to H(k) =n × Smatch − k × Ge when processing k operators. The
parameter p is given such that H(p) ≥ 0, which results in
Eq. (11).
We implemented this idea by extending a GPUaccelerated string search implementation [9] that localized
scan computation. Figure 4 shows how GPU threads process
local scan in parallel. The local scan is processed in four
steps. Without loss of generality, 32-bit data here is assumed
to simplify the description.
1)

2)

3)

4)

Scan within segments. For example, the third
thread block in Fig. 4 computes a vector(x8 , x8-9 , x8-10 , x8-11 ) for its responsible segment(x8 , x9 , x10 , x11 ).
Gather representatives from segments. To perform
scan among segments, we collect the last elementsfrom locally scanned segments. In Fig. 4, four representatives x0-3 , x4-7 , x8-11 , and x12-15 are copied
to a temporal buffer.
Scan for gathered values. Similar to the firststep, local scan is applied to the temporalbuffer. After that, the buffer holds scanned results(x0-3 , x0-7 , x0-11 , x0-15 ), each starting from the initial element. This step can be omitted when thesegment size is greater than p.
Scatter scanned values to segments. The scanned
results obtained at the third step are scattered toall segments. The scattered results are then added
to the local segments. Finally, every element holds
the scanned result starting from the initial element.

x0

x1

x2 x3

▷ initialize with 8-bit mode
▷ for each row▷ update the highest score

x4 x5 x6 x7

x8 x9 x10 x11 x12 x13 x14 x15

1. Scan within segments
x0 x0-1 x0-2 x0-3 x4 x4-5 x4-6 x4-7 x8 x8-9 x8-10 x8-11 x12 x12-13 x12-14 x12-15

2. Gather representatives
from segments

x0-3

x4-7

x8-11 x12-15

x0-3

x0-7

x0-11 x0-15

ӻӻ

3. Scan for gathered valuesӻ

4. Scatter scanned values
to segments
x0 x0-1 x0-2 x0-3 x4 x4-5 x4-6 x4-7 x8 x8-9 x8-10 x8-11 x12 x12-13 x12-14 x12-15

x0 x0-1 x0-2 x0-3 x0-4 x0-5 x0-6 x0-7 x0-8 x0-9 x0-10 x0-11 x0-12 x0-13 x0-14 x0-15

Figure 4. Localized parallel scan on the GPU. The input array is scannedin four steps. Barrier synchronization is required between steps.

4. Experimental Results
We evaluated the proposed method in terms of execution time and compared the method with previous GPUaccelerated methods [3], [13]. We measured the alignment
throughput in cell updates per second (CUPS), given by
mn/t, where t is the execution time that contains (1) thetime t1 spent on the GPU, (2) that t2 on the CPU, and (3)that t3 needed for data transfer between the main memory
and the GPU memory.
For experiments, we used a desk-top PC equipped withan Intel Core i7-6700 3.4 GHz CPU, 32 GB main memory,and an NVIDIA GeForce GTX 980ti GPU with 6 GB device
memory. All implementations ran on Windows 10. The

TABLE 1. B REAKDOWN OF EXECUTION TIME FOR PROPOSED AND
PREVIOUS METHODS [13] WITH n = 256 . t1 , t2 , AND t3 REPRESENTCPU TIME , GPU TIME , AND CPU-GPU TRANSFER TIME INMILLISECONDS , RESPECTIVELY.

Throughput (GCUPS)

12

Our method32-bit scan [13]10Anti-diagonal [3]Breakdown

8

t1
t2
t3
Total

6
4
2
0
1,024

8,192

65,536524,288

4,194,304

33,554,432

Throughput (GCUPS)(a)
14

Throughput (GCUPS)

Our method32-bit scan [13]10Long (m = 33M)
Proposed Previous
670.0
2217.0
7.1
6.2
30.0
2.3
707.1
2225.500 00 00 00 00 00 00 00 00 00 00
000000000000
Our method 0 0 0 0 0 0 0 0 0 0 0 0 32-bit scan [13]m : Length of long sequence

12

Short (m = 100K)
Proposed Previous
37.0
44.0
3.3
1.2
13.0
0.3
53.3
45.5

Anti-diagonal [3]

8

20
18
16
14
12
10
8
6
4
2
0

100
80

00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00 00 00 00 00
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 0 0 0 0 0 0
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
000000200

6Ratio of 8-bit mode

400

600

800

60
40
20

r : Ratio of 8-bit mode (%)

14

0

1000

n : Length of short sequence

4

Figure 6. Alignment throughput with different ratio r of 8-bit mode over
32-bit mode. r = 0 and r = 100 mean that the alignment is fully processedwith 32-bit mode and 8-bit mode, respectively.

2
0
1,024

8,192

65,536524,288

4,194,304

33,554,432m : Length of long sequence(b)
Figure 5. Alignment throughput with different sequence length m. Resultsfor (a) n = 128 and (b) n = 256. The horizontal axis is given in alogarithmic scale.

methods were implemented using Visual Studio 2013 andCUDA 7.5 [10].
As for long sequences, we used the whole genome
shotgun sequence of Homo sapiens chromosome 21 (accession number: [GenBank CM000511]). For short sequences,
we used a segment of Influenza A virus gene (accession
number: [GenBank FJ966086]). Alignments were carriedout with a affine gap penalty (Gs = 5 and Ge = 1) and ascoring matrix S such that Smatch = 1 and Sunmatch = −2.
Figure 5 shows the measured alignment throughputs
with different lengths m of sequences ranging from 1 Kbp
to 33 Mbp. Our method was faster than the previous 32bit scan-based method [13] when m > 524K. For such
short sequences, i.e., short vs short sequences, the overheadof data transfer limited the alignment throughput of our
method. However, the speedup over the previous method
reached a factor of 3.1 when n = 256 and m = 33M. Furthermore, the proposed method was 21.2 times faster than
the antidiagonal-based method [3]. The highest throughput
of 12.3 GCUPS was obtained when n = 256 and m = 33M.

We next investigated the breakdown of execution time t
for short and long sequences (Table 1). When m = 100K
and n = 256, i.e., for short sequences, our method was
slower than the scan-based previous method. In contrast,
our method was faster than the previous method for longsequences. For both cases, our method achieved shorter GPU
time t1 . However, the CPU-GPU transfer time t3 increased
from 2.3 ms to 30.0 ms when m = 33M and from 0.3 ms
to 13.0 ms when m = 100K. This increase was due to thedata transfer occurred at line 6 in Algorithm 1. That is, theproposed method transferred the highest scores to switch thedata structure according to Eqs. (9) and (10) computed on
the CPU.
Finally, we counted the number l of matrix elements
processed with the 8-bit mode. According to this numberl, we computed r = l/mn, i.e., the rate of the 8-bit modeover the 32-bit mode (Fig. 6). To investigate the impactof the 8-bit mode, we varied the rate r by aligning the
genome shotgun sequence of Homo sapiens chromosome
21 with itself; we used the sequence as a long sequence
while its subsequence as a short sequence of 200–1000 bp.
In Fig. 6, the proposed method shows a linear relationshipbetween the ratio r and the throughput. That is, the achieved
acceleration over the previous method comes from the 8-bit
data structure, which reduces the amount of memory access.
In particular, the proposed method is effective when n issmall or the sequences are not highly similar.

5. Conclusion
We have presented a mixed-word-length method capableof accelerating the SW algorithm for a pair of short andlong sequences. The key idea for realizing acceleration
is to integrate an 8-bit data structure into a scan-based
parallel scheme [13]. Our method not only saves the memory
bandwidth but also avoids integer overflow by appropriately
switching 8-bit mode and 32-bit mode.
Experimental results show that the proposed method isthree times and twenty times faster than the scan-based
method and the antidiagonal-based method, respectively.
Therefore, the proposed method is useful for accelerating
SW alignment for a pair of short and long sequences.
Future work includes further optimization such as pruning of scoring computation and reduction of memory access.

ACKNOWLEDGMENTS
This study was supported in part by the Japan Society
for the Promotion of Science KAKENHI Grant Numbers
15H01687, 16H02801 and 15K12008.

References[1]T. F. Smith and M. S. Waterman, “Identification of common molecularsubsequences,” J. Molecular Biology, vol. 147, no. 1, pp. 195–197,
Mar. 1981.

[2]

F. Ino, Y. Munekawa, and K. Hagihara, “Sequence homology searchusing fine grained cycle sharing of idle GPUs,” IEEE Trans. Paralleland Distributed Systems, vol. 23, no. 4, pp. 751–759, Apr. 2012.

[3]

M. Korpar and M. Šikić, “SW#–GPU-enabled exact alignments on
genome scale,” Bioinformatics, vol. 29, no. 19, pp. 2494–2495, Oct.
2013. [Online]. Available: https://sourceforge.net/swsharp/code/[4]

E. F. de O. Sandes, G. Miranda, A. C. M. A. de Melo, X. Martorell, and E. Ayguadé, “CUDAlign 3.0: Parallel biological sequencecomparison in large GPU clusters,” in Proc. 14th IEEE/ACM Int’l
Symp. Cluster, Cloud and Grid Computing (CCGrid’14), May 2014,
pp. 160–169.

[5]

D. Okada, F. Ino, and K. Hagihara, “Accelerating the Smith-Waterman
algorithm with an interpair pruning method for all-pairs comparisonof base sequences,” BMC Bioinformatics, vol. 16, no. 321, Oct. 2015,15 pages.

[6]

E. Rucci, C. Garcia, G. Botella, A. D. Giusti, M. R. Naiouf, and
M. Prieto-Matias, “First experiences accelerating Smith-Waterman
on Intel’s Knights Landing processor,” in Proc. 17th Int’l Conf.
Algorithms and Architectures for Parallel Processing (ICA3PP’17),
Aug. 2017, pp. 569–579.

[7]

J. D. Owens, M. Houston, D. Luebke, S. Green, J. E. Stone, and J. C.
Phillips, “GPU computing,” Proceedings of the IEEE, vol. 96, no. 5,
pp. 879–899, May 2008.

[8]

K. Ikeda, F. Ino, and K. Hagihara, “Efficient acceleration of mutual information computation for nonrigid registration using CUDA,” IEEE
J. Biomedical and Health Informatics, vol. 18, no. 3, pp. 956–968,
May 2014.

[9]

Y. Mitani, F. Ino, and K. Hagihara, “Parallelizing exact and approximate string matching via inclusive scan on a GPU,” IEEE Trans.
Parallel and Distributed Systems, vol. 28, no. 7, pp. 1989–2002, Jul.
2017.

[10] NVIDIA Corporation, “CUDA C Programming Guide Version
8.0,” Jan. 2017. [Online]. Available: http://docs.nvidia.com/cuda/pdf/
CUDA C Programming Guide.pdf[11] W. Liu, B. Schmidt, G. Voss, and W. Müller-Wittig, “Bio-sequence
database scanning on a GPU,” in Proc. 20th IEEE Int’l Parallel andDistributed Processing Symp. (IPDPS’06), Apr. 2006.
[12] Y. Munekawa, F. Ino, and K. Hagihara, “Design and implementationof the Smith-Waterman algorithm on the CUDA-compatible GPU,”in Proc. 8th IEEE Int’l Conf. Bioinformatics and Bioengineering(BIBE’08), Oct. 2008, 6 pages (CD-ROM).
[13] A. Khajeh-Saeed, S. Poole, and B. Perot, “Acceleration of the SmithWaterman algorithm using single and multiple graphics processors,”
J. Computational Physics, vol. 229, no. 11, pp. 4247–4258, Jun. 2010.
[14] W. Liu, B. Schmidt, G. Voss, and W. Müller-Wittig, “Streaming
algorithms for biological sequence alignment on GPUs,” IEEE Trans.
Parallel and Distributed Systems, vol. 18, no. 9, pp. 1270–1281, Sep.
2007.
[15] N. Bell and J. Hoberock, Thrust: A Productivity-Oriented Library
for CUDA. San Mateo, CA: Morgan Kaufmann, Jan. 2011, ch. 26,
http://thrust.github.io/.