Parallel Computing 116 (2023) 103019 Contents lists available at ScienceDirect Parallel Computing journal homepage: www.elsevier.com/locate/parco GPU acceleration of Levenshtein distance computation between long strings David Castells-Rufas โˆ— Department of Microelectronics & Electronic Systems, Universitat Autรฒnoma de Barcelona, Cerdanyola del Vallรฉs, 08193, Spain Barcelona Supercomputing Center, Barcelona, 08034, Spain ARTICLE INFO ABSTRACT Keywords: Computing edit distance for very long strings has been hampered by quadratic time complexity with respect Levenshtein distance to string length. The WFA algorithm reduces the time complexity to a quadratic factor with respect to the Edit distance edit distance between the strings. This work presents a GPU implementation of the WFA algorithm and a WFA algorithm new optimization that can halve the elements to be computed, providing additional performance gains. The GPU implementation allows to address the computation of the edit distance between strings having hundreds of Parallel processing millions of characters. The performance of the algorithm depends on the similarity between the strings. For strings longer than million characters, the performance is the best ever reported, which is above TCUPS for strings with similarities greater than 70% and above one hundred TCUPS for 99.9% similarity. 1. Introduction using the WFA algorithm with O (๐‘‘) space complexity on Graphical Processing Units (GPUs) for fast execution times and limited memory The Levenshtein distance [1] measures the number of edit op- use. erations (insertion, deletion, and substitution) required to convert a This paper is organized as follows. In Section 2 we analyze previous string P into a string T. Levenshtein distance is also known as edit work. In Section 3 we analyze the fundamentals of the WFA algorithm distance (ED). It is used in several computing domains, including text and its adaptation to the computation of Levenshtein distance. The analysis [2], music analysis [3], and bioinformatics [4]. Computer adaptation to GPU architectures is presented in Section 4, where we Science uses the concepts of time complexity and space complexity [5] follow an iterative improvement process to get the maximum perfor- to be able to compare algorithms. These concepts refer to the relation mance. Finally, we analyze the results and contrast them with other between the problem size (๐‘›) and the execution time and the required works from the state of the art in Section 5 before concluding. memory, respectively. Wagnerโ€“Fischer (WF) [6] proposed a dynamic programming (DP) algorithm to compute ED with an O (๐‘š ร— ๐‘›) time 2. Previous work and space complexity, where ๐‘š = |๐‘ƒ | and ๐‘› = |๐‘‡ | (i.e. the lengths of the strings). The space complexity is given by the need to store the DP The classic WF algorithm is based on building a DP table ๐ท using table to be able to back-trace the alignment sequence between the two strings. If only the distance score is required, time complexity remains Eq. (1). the same but space complexity can be reduced to O (๐‘š). โŽง๐‘– , if ๐‘— = 0 The Wavefront Alignment Algorithm (WFA) is a DP technique that โŽช ( ) โŽช๐‘— , if ๐‘– = 0 reduces the time and space complexity to O ๐‘‘ 2 , where ๐‘‘ is the โŽช edit distance between both strings. WFA was originally proposed to ๐ท๐‘–,๐‘— =โŽจ โŽง๐ท๐‘–โˆ’1,๐‘— + 1 (1) โŽชmin โŽช , otherwise compute Smith-Waterman (SW) gap-affine global alignment [7] but it โŽช โŽจ๐ท๐‘–,๐‘—โˆ’1 + 1 can be adapted to compute other distance or similarity scores based on โŽช โŽช โŽฉ โŽฉ๐ท๐‘–โˆ’1,๐‘—โˆ’1 + [๐‘ƒ [๐‘– โˆ’ 1] โ‰  ๐‘‡ [๐‘— โˆ’ 1]] computing a DP table. In [8], it was adapted to compute Levenshtein distance for limited strings lengths up to a thousand characters. The The ED between ๐‘ƒ and ๐‘‡ is determined by the cell ๐ท๐‘š,๐‘› from the problem of using DP algorithms on very long strings has historically table. This algorithm can be considered a brute-force approach where been hindered by the difficulty of managing the high memory demand all possible edit operations are considered but it uses memoization to and excessive computing time. This paper proposes several optimiza- avoid repeating the same computing for several results. Eq. (1) can be tions to compute the Levenshtein distance between very long strings found in the literature in different ways. In our case, we use the Iverson โˆ— Correspondence to: Department of Microelectronics & Electronic Systems, Universitat Autรฒnoma de Barcelona, Cerdanyola del Vallรฉs, 08193, Spain. E-mail address: david.castells@uab.cat. https://doi.org/10.1016/j.parco.2023.103019 Received 4 October 2022; Received in revised form 16 March 2023; Accepted 30 March 2023 Available online 3 April 2023 0167-8191/ยฉ 2023 The Author. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). D. Castells-Rufas Parallel Computing 116 (2023) 103019 bracket notation for the extra value when characters do not match. In and multi-threading opportunities available on modern CPUs to get Iverson bracket notation [9], if the expression inside the brackets is true performance above 100 GCUPS. the value is considered 1, and zero otherwise. On the other hand, GPUs offer a massively parallel execution plat- The ED is not the only interesting result that can be obtained form. DP algorithms have been adapted to GPUs in search of scalability. with the WF algorithm. Many applications recover the sequence of In [15], they reached a significant performance speed (above 1 TCUPS) operations required to transform P into ๐‘‡ by analyzing the elements using the Myers BPM algorithm on strings up to 1k characters. Longer of the ๐ท table. This sequence can be represented as an alignment sequences (8k and 22.5k) were targeted by Balhalf in [16,17] using an string, where M is used to denote a match, and I,D, and X are used anti-diagonal parallel approach with modest throughput results under 1 to denote insert, delete, and substitution operations respectively. There GCUPS. Tiling [18] is a general strategy to increase data locality in GPU can be multiple valid alignment strings, as different edit operations workloads. It is used in [19,20], and [21] to address longer strings of up can result into the same ED. The following example illustrates one to 2 M with a peak performance of 50 GCUPS. FPGAs also offer a good alignment string obtained from two short strings, implementation platform for applications with a regular data layout. Our recent work [22] achieves more than 5 TCUPS for short genomics strings (<300 characters). Computer clusters have not been generally P: -ACCATGGACTG- used to speed up this problem. Sadiq and Yousaf [23] is one of the few ||| || ||| works trying to parallelize ED for clusters with modest results. T: CACC-TG-ACTTA The low research attention given to ED between long strings con- alignment: IMMMDMMDMMMXI trast with the enormous interest of other related DP algorithms such distance: 5 as SW [24] and their variants. The reason for this is its superlative impact in bioinformatics, especially genomics, which always demands faster methods to compare long DNA strings. The advances in this The backtrace method to obtain the alignment string is based on algorithm are making possible genome-wide string comparisons. The recovering the active branches of the conditional statements of Eq. (1). human genome length is approximately 3 G. Many parallel strategies When dealing with very long strings, the cost of storing the table is have been studied to compare at that scale [25]. prohibitive, and piece-based strategies (such as the one proposed by Some of the most relevant results on accelerating SW for long Hirschberg [10]) are required to be able to obtain an alignment string. strings with GPUs have been achieved by the several versions of the This work focuses on obtaining the distance between the input strings. CUDAlign library [26โ€“29]. Their goal is to recover the full alignment. Obtaining the alignment string is left for future works. Their approach is based on sampling some of the rows during DP table There have been many attempts to improve the performance of computation to later refine the alignment in a piece-wise approach matrix-based DP algorithms. The DP matrix offers a very regular layout using massive parallelism. A large part of the execution time is devoted that seems adequate for a parallel implementation at a first glance. to the first phase of the algorithm, which is the building of the DP However, the conditional branches appearing in Eq. (1) imply an table. Its complexity is comparable to computing the ED. Thus, we will irregular data dependency among the cells, and makes the contribution use datasets from SW works to compare the execution time of our ED of many of the computed cells irrelevant to the final result. Hence, solution. the methods to speed up the algorithm are focused on two aspects: (1) speeding up the process computation of the cells. (2) avoiding 3. WFA algorithm computing some of the cells of the DP matrix. In the first group, parallel processing is usually applied. The ex- The wavefront algorithm (WFA) [7] belongs to the second group of isting data dependency requires to use of an anti-diagonal approach algorithms, which avoid to compute cells to speed-up the computation which is hard to parallelize since multiple threads receive an irregular of the DP table. WFA only computes the cells from the diagonals of the workload. On the other hand, other approaches target the speedup DP table where the distance is incremented. As the algorithm is based of the single-threaded execution. The BPM Myers algorithm [11] is a on diagonals, it is convenient to define a coordinate system that refers remarkable example. It uses bit-vectors and lookup tables to reduce to the elements of the matrix ๐ท in terms of diagonals and displacements the computational load by a factor related to the word length of the within the diagonals. computing architecture. In the second group, Ukkonen [12] proposed Several coordinate systems could be used to identify the cells of the a banded approach, which was later combined with Myersโ€™ bit vectors matrix as a displacement within the diagonals. In Fig. 1 we show the by Hyyrรถ [13]. four coordinate systems we foresee: offset, radius, row, and column. We These two different strategies make comparing the performance of can describe the ๐ท table using any of these coordinate systems. In the different works tricky. We know that DP algorithmsโ€™ time complexity is original Marco-Solaโ€™s paper [7] a column coordinate system is used. In O(๐‘š ร— ๐‘›). Hence, many works report the number of computed cells from the rest of this paper we will use a radius coordinate system, which the DP table in cell updates per second (CUPS), which are defined as might be more intuitively understood by the reader. Using a different Eq. (2). coordinate system has an impact on the definition of the equations used ๐‘šร—๐‘› CUPS = (2) by the algorithm, which we will introduce below. However, they can ๐‘ก be easily adapted to all coordinate systems. However, when using algorithms from the second group, the CUPS The transformation between matrix (Cartesian) ๐ท and pyramid unit could be misleading as they do not require computing all cells but (radius) ๐ทฬ‚ coordinate systems is described with Eqs. (3) and (4), which a subset of them. Nevertheless, we think that CUPS are still valid as a also use the Iverson bracket notation. An example ๐ท table with its throughput measure as it maintains an inverse relation (i.e. 1โˆ•๐‘ฅ) with corresponding ๐ท ฬ‚ pyramid are shown in Figs. 2(a) 2(b) respectively. execution time. The Edlib library [4] is a popular implementation of the BPM ฬ‚๐‘—โˆ’๐‘–,[๐‘–>๐‘—]๐‘—+[๐‘–โ‰ค๐‘—]๐‘– ๐ท๐‘–,๐‘— = ๐ท (3) algorithm, which exploits the intrinsic parallelism of the add instruction in carry propagation to compute several cells from the table in a single ฬ‚๐‘˜,๐‘Ÿ = ๐ท[๐‘˜<0](โˆ’๐‘˜)+๐‘Ÿ,[๐‘˜โ‰ฅ0](๐‘˜)+๐‘Ÿ ๐ท (4) instruction execution. With the end of Denardโ€™s scaling regime, comput- ing platforms have been forced to embrace parallel architectures. The We define the wavefront pyramid ๐‘Š ฬ‚ as the pyramid that contains BGSA [14] is a remarkable extension to benefit from the vectorization ฬ‚ in ๐‘Š๐‘˜,๐‘Ÿ the farthest displacement in the diagonal ๐‘˜ of the original 2 D. Castells-Rufas Parallel Computing 116 (2023) 103019 3.1. Algorithm optimizations derived from data dependency In the WF algorithm, the solution to the ED problem is found in ๐ท๐‘š,๐‘› . However, in WFA, we can only ensure that the solution is in the diagonal ๐‘˜๐‘ ๐‘œ๐‘™ = ๐‘›โˆ’๐‘š, and that its radius is ๐‘Ÿ๐‘ ๐‘œ๐‘™ โ‰ค ๐‘š๐‘Ž๐‘ฅ(๐‘š, ๐‘›). From Eq. (7), it is obvious that ๐ถฬ‚๐‘˜๐‘ ๐‘œ๐‘™ ,๐‘Ÿ only depends on the neighbors of the cell in previous radius (i.e. ๐ถฬ‚๐‘˜ โˆ’1,๐‘Ÿโˆ’1 , ๐ถฬ‚๐‘˜ ,๐‘Ÿโˆ’1 , and ๐ถฬ‚๐‘˜ +1,๐‘Ÿโˆ’1 ). These cells, in ๐‘ ๐‘œ๐‘™ ๐‘ ๐‘œ๐‘™ ๐‘ ๐‘œ๐‘™ turn, will depend on their neighbors from the previous radius, and so forth. Therefore, in the worst case, the only required cells to compute the ๐ถฬ‚ pyramid, and hence, the ๐‘Šฬ‚ pyramid, will be the cells inscribed inside the diamond shaped parallelogram described by Eq. (9). ๐‘Ÿ โˆˆ [0, ๐‘Ÿ๐‘ ๐‘œ๐‘™0 ] (9) ๐‘˜ โˆˆ [๐‘š๐‘–๐‘›(๐‘Ÿ, ๐‘˜๐‘ ๐‘œ๐‘™0 โˆ’ (๐‘Ÿ๐‘ ๐‘œ๐‘™0 โˆ’ ๐‘Ÿ)), ๐‘š๐‘–๐‘›(๐‘Ÿ, ๐‘˜๐‘ ๐‘œ๐‘™ + (๐‘Ÿ๐‘ ๐‘œ๐‘™0 โˆ’ ๐‘Ÿ))] Unkkonnen and later Myers [11] and Hyyrรถ [13] exploited the fact that, in ED, the increments on the diagonals can only be 0 or +1. In the worst case (having no match between ๐‘ƒ and ๐‘‡ ), the values of the solution diagonal in both ๐ท ฬ‚ and ๐‘Šฬ‚ will increase from 0 to ๐‘Ÿ๐‘ ๐‘œ๐‘™ in +1 increments. In case of a single match in the diagonal, all increments will be +1 except a single value being zero, and then, the resulting ED will be ๐‘Ÿ๐‘ ๐‘œ๐‘™ โˆ’ 1. After this observation, we can generally state that after computing a value of the solution diagonal, ๐‘“ = ๐‘Š ฬ‚๐‘˜ ,๐‘Ÿ , the new worst ๐‘ ๐‘œ๐‘™ case radius must be ๐‘Ÿ๐‘ ๐‘œ๐‘™๐‘Ÿ โ‰ค ๐‘Ÿ๐‘ ๐‘œ๐‘™0 โˆ’ (๐‘“ โˆ’ ๐‘Ÿ). We can use the new upper bound to dynamically define the range of cells to compute from the ๐‘Š ฬ‚ pyramid as the algorithm progresses along the pyramid columns. This approach, that we name dynamic diamond, reduces even more the number of the required cells to compute. Fig. 1. Possible diagonal coordinate systems using an index to identify the diagonal In the example illustrated in Fig. 2, 36 cells from ๐‘Š ฬ‚ are computed and another index to identify the displacement within the diagonal. In the figures, the arrows define the diagonal coordinates, and the values on the boxes refer to the to get the final result using the classic WFA algorithm illustrated by displacement value within the diagonal. Fig. c. With the Dynamic Diamond approach illustrated by Fig. d, we require to compute 24 cells, which represents a factor 1.5 ร— of savings. Actually, the savings of this approach will be between 1ร— and 2ร—, being ฬ‚ pyramid having ED ๐‘Ÿ. ๐‘Š ๐ท ฬ‚ can be formally derived from ๐ท ฬ‚ using higher as the distance between ๐‘ƒ and ๐‘‡ is higher. If the alignment path is not required, we can reduce the storage Eq. (5). to just two columns from the ๐‘Š ฬ‚ pyramid, as the whole pyramid is ฬ‚๐‘˜,๐‘Ÿ = max ๐‘ž โˆฃ ๐ท ๐‘Š ฬ‚๐‘˜,๐‘ž = ๐‘Ÿ (5) only required to perform the backtrace. In this case, when computing โˆ€๐‘ž a column, only the information from the previous column as can be The great advantage of the WFA algorithm is that it finds an alter- observed by the data dependency pattern. We will have to dimension ฬ‚ pyramid without computing ๐ท. native way of building the ๐‘Š ฬ‚ Instead, the column as the worst case, which corresponds to the final column of it uses two additional pyramids ๐ถฬ‚ and ๐ธ ฬ‚ as described in Eq. (6). These the pyramid containing ๐‘š + ๐‘› + 1 elements. pyramids are described below. 4. GPU implementation ฬ‚๐‘˜,๐‘Ÿ = ๐ถฬ‚๐‘˜,๐‘Ÿ + ๐ธ ๐‘Š ฬ‚ ฬ‚ (6) ๐‘˜,๐ถ๐‘˜,๐‘Ÿ GPUs offer a massively parallel execution platform based on a The compute pyramid ๐ถฬ‚ is built by selecting the maximum values large number of processing elements (PEs) and a multi-ported memory ฬ‚ pyramid as described in Eq. (7). from the depending cells of the ๐‘Š hierarchies with high memory bandwidth. One difference of GPU archi- tectures with respect to general purpose chip multiprocessors (CMP) is โŽง0 , if ๐‘˜ = 0 โˆง ๐‘Ÿ = 0 โŽช that, at a given time, a large number (if not all) PEs execute the same โŽช โŽง๐‘Š๐‘˜+1,๐‘Ÿโˆ’1 + [๐‘˜ < 0] program (called kernel), which is typically small in size. By doing this, ๐ถ๐‘˜,๐‘Ÿ =โŽจ โŽช (7) โŽชmax โŽจ๐‘Š๐‘˜,๐‘Ÿโˆ’1 + 1 , otherwise instruction fetch can be shared among a large number of PEs in a SIMD โŽช โŽช fashion. To maximize the efficiency of instruction execution among โŽฉ โŽฉ๐‘Š๐‘˜โˆ’1,๐‘Ÿโˆ’1 + [๐‘˜ > 0] the PEs, control flow divergence is minimized with several techniques, The extend pyramid ๐ธ ฬ‚ is the pyramidal form of the extend table ๐ธ such as predicated conditionals, stack memory avoidance, low latency as defined in Eq. (8), which contains the length of the longest common context switching, and large number of in-flight threads and memory prefixes of P and ๐‘‡ starting at positions ๐‘– and ๐‘— respectively. access latency hiding techniques. Several programming frameworks can be used to program GPUs. ๐ธ๐‘–,๐‘— = ๐ฟ๐ถ๐‘ƒ (๐‘ƒ๐‘– , ๐‘‡๐‘— ) (8) NVIDIA favors CUDA and OpenACC but it also supports OpenCL, which offers a wider platform support for other accelerator technologies like ฬ‚ are usually needed In practice, only a small subset of values from ๐ธ FPGAs. Our implementation is mainly based on OpenCL, although as the radius of the accessed elements depend on the values of the we also code some variants in CUDA because of the richer available previously computed elements of ๐ถฬ‚ pyramid as shown in Eq. (6). profiling tools. Before going further into the implementation, we should Moreover, the data dependency allows to avoid storing the pyramids introduce some GPU programming concepts. ฬ‚ ๐ธ, ๐ถ, ฬ‚ which are computed as needed. Fig. 2 depicts an example of a ๐ท The GPU has a global memory external to the chip. This is usually a matrix, its pyramidal form ๐ท,ฬ‚ and the wavefront pyramid ๐‘Š ฬ‚. multi-ported big memory implemented in several banks of some of the 3 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Fig. 2. Visual example of how the classic DP approach is transformed into the WFA approach and the dynamic diamond optimization presented in this paper. (a) DP 12 ร— 12 table ๐ท to compute the distance between the example strings P = โ€˜โ€˜CACCTGACTTAโ€™โ€™ and T = โ€˜โ€˜ACCATGGACTGโ€™โ€™. (b) Pyramidal form ๐ท ฬ‚ of the DP table in radius coordinates. (c) Pyramidal form of the WFA algorithm. (d) Dynamic-diamond optimization presented in 4.3 to reduce the number of computed cells. The arrows represent the farthest pyramid column that the algorithm considers. As the algorithm computes more columns from left to right, the furthest pyramid column is reduced and less cells must be computed thanks to the diamond shaped data dependency. Table 1 Equivalence among different terminology in OpenCL and CUDA. Concept OpenCL CUDA Processor Compute unit Streaming multiprocessor Workload unit Workitem Thread Group of workload units executed in the same processor Workgroup Thread block Batch of workload units computed in a kernel invocation NDRange Grid latest DRAM technologies. โ€˜โ€˜Single Instruction Multiple Dataโ€™โ€™ (SIMD) processing elements (PEs) are grouped into computing units (CU). NVIDIA calls them Streaming Multi-processors (SM). PEs in the same CU can be synchronized by means of specific synchronization primitives and they can share some on-chip memory. When invoking kernels in a GPU device, the host specifies the number of workitems required to execute the workload. A workitem is a unit of work associated to an invocation of the kernel with a unique identifier that can be a 1D, 2D, or 3D index. Workgroups are groups of workitems, and they can be used to ensure that a range of workitems are executed in the same CU. A scheduler distributes the workitems among the available PEs. Once a workitem starts executing in a PE it can be context-switched, but it generally will not be retired from the CU until it has completed its Fig. 3. Data dependency to compute a cell (a) from the ๐‘Š ฬ‚ pyramid. For each cell, execution. Contrary to what happens in a normal preemptive multitask- three cells from the previous column are accessed (b,c,d), as expressed in Eq. (7). One ing system, there might be workitems (threads) not starting until there or more characters from the ๐‘ƒ and ๐‘‡ strings are also accessed to compute the LCP (corresponding to Eq. (8)). are available PEs for execution. This behavior prevents from having global synchronization primitives among all workitems participating in a kernel invocation. Global synchronization can only occur by host kernel invocations, as the host can ensure that the whole batch of by Eq. (10), where ๐›ด is the alphabet of the strings being compared. workitems launched are completed and cache memories are flushed. ( )๐‘ข 1 The CUDA framework from NVIDIA uses a different terminology, but ๐‘ƒ (run ๐‘ข) โ‰ˆ (10) |๐›ด| the equivalence is shown in Table 1. As stated in the previous section, WFA space complexity can be Thus, in practice, the majority of accesses will require very few iter- reduced by just storing two columns of the pyramid: the currently being ations during the extension phase, which will result in a low divergence computed and the previous one. An easy GPU parallelization strategy is between workitems. Moreover, we expect a large number of cache hits to devote one workitem to every cell of the pyramid column currently as workitems access nearby positions. computed. As depicted in Fig. 3, each workitem computing a cell must When working with short strings, the workload is relatively small access three cells values (๐‘, ๐‘, ๐‘‘) from the previous column to obtain compared to the large number of processing elements (PEs) available on the ๐ถฬ‚๐‘‘,๐‘Ÿ value and must access one or several values from the ๐‘ƒ and ๐‘‡ a GPU, resulting in low GPU occupancy. To address this, it is common strings to finally assign the ๐‘Ž value to the ๐‘Šฬ‚๐‘‘,๐‘Ÿ pyramid. practice to group a large number of short string pairs and process The probability to have a run of ๐‘ข consecutive matches in the sub- them in parallel, as described in [8]. However, when dealing with strings from ๐‘ƒ and ๐‘‡ during the extension phase can be approximated large strings, even with a pyramidal data structure, the amount of work 4 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Fig. 4. Performance achieved in NVIDIA GeForce RTX 3090 with different kernel invocation strategies as a function of the workgroup size. GCUPS are computed as defined in Eq. (2). (a) Baseline. (b) Adapting the number of workitems as the algorithm progresses along the columns has a beneficial impact in performance. (c) However, the biggest additional benefit is achieved when we reduce the frequency at which the host reads the memory positions that indicate the processing is complete. required to compute a single column can far exceed the available PEs is computed to start the next one. We use kernel invocation as a global after the initial iterations. synchronization method. For every column, we invoke the kernel with Several strategies can be used trying to maximize the throughput an NDRange to compute all the possible cells of the column, which of the system. In the following subsections we will iteratively analyze correspond to the worst case (๐‘š + ๐‘› + 1 elements). Each workitem some implementation alternatives and measure their impact on the determines whether the cell has to be executed or not. In the later case, system performance. To demonstrate these approaches, we will use two it finishes early. files from the National Center for Biotechnology Information (NCBI), Although there is no local memory sharing or fine-grain synchro- with each file being over 3 million characters long, corresponding to nization among workitems, the workgroup size has an influence in the accessions BA000035.2 and BX927147.1. The ED between these files execution time, since the GPU benefits from having multiple workitems exceeds 1.5 million, providing an opportunity to achieve significant executing the same instruction in a SIMD fashion and workitems are performance improvements. We will report performance in GCUPS. scheduled in multiples of the SIMD lanes of the GPU. Fig. 4(a) shows the achieved performance as a function of the workgroup size. The 4.1. Kernel invocation strategies workgroup size is known to influence performance and exhibit a zig- zag pattern [30]. The fluctuations are explained by changes in the Memory requirements are some of the fundamental aspects that amount of used hardware and occupancy of the SIMD units. The influence the performance of an implementation. In the worst case sce- figure considers the total execution time of the application, not just nario of the strings being completely different (i.e. not a single match), a single kernel invocation. The maximum performance (156 GCUPS) is the size of the required pyramid is the same than the DP table, i.e. ๐‘›2 achieved with a workgroup size of 380 workitems, and corresponds to for same-length strings. Therefore, the supported maximum โˆš length of an execution time of 65.8 s. In this strategy, the host invokes the kernel same-length strings is determined by the expression ๐‘›๐‘š๐‘Ž๐‘ฅ < ๐‘€๐บ๐‘ƒ ๐‘ˆ โˆ•๐‘ค, as many times as the number of required columns from the ๐‘Š ฬ‚ pyramid. where ๐‘€๐บ๐‘ƒ ๐‘ˆ is the global memory size of the GPU and ๐‘ค is the size We analyze how the kernel execution varies as the algorithm pro- of the data type used to store the pyramid values (both measured in gresses in column radius. The blue line in Fig. 5 depicts the execution bytes). In our case, we use an NVIDIA RTX 3090 with ๐‘€๐‘…๐‘‡ ๐‘‹3090 = 24 time of kernel invocations. Although invocations are launched with the GB and a ๐‘ค = 4 B, so ๐‘›๐‘š๐‘Ž๐‘ฅโˆ’๐‘…๐‘‡ ๐‘‹3090 < 77 k. Therefore, using the standard worst case number of workitems (the rows from the pyramid column), WFA algorithm, the pyramid for longer strings would not fit into the early detection of useless computation in the OpenCL kernel code will GPU memory. For longer strings, we have to use the two column make many threads to exit in few instructions. The ๐‘ฅ axis of the figure computation optimization presented at the end of Section 3.1. In the represents the radius of the column in the kernel invocation. For low general case, we are not restricted to same-length strings. We allocate values, the execution time remains constant and many workitems exit space in global memory to store two columns of length ๐‘š + ๐‘› + 1, which by checking early exit conditionals. At a certain threshold the number corresponds to the maximum possible number of rows of the pyramid. of useful workitems are involved in slightly more complex code that Note that we generally use a fraction of the cells of the two columns, increases the overall execution time. As the cells invocations reach the since the full column will be completely used in the right-most column end point many workitems from distant diagonals will already reach of the pyramid. Although dynamically allocating memory required by early exit conditions making a decrease on the overall execution time each column as the algorithm progresses may be considered as a space- of those kernel invocations. This approach obviously has an overhead, saving measure, it is counterproductive due to the requirement for the since it launches more workitems than are really necessary. host to coordinate the memory allocation process, which is known to be slow. With the fixed size strategy, the number of cells required for 4.1.2. Dynamic number of workitems same-length strings is 2(2๐‘›+1). So, in this case the maximum supported In a second implementation we invoke the kernel with a dynamic string length would be determined by the expression ๐‘›๐‘š๐‘Ž๐‘ฅ < ๐‘€๐บ๐‘ƒ ๐‘ˆ โˆ•4๐‘ค. number of workitems in the NDRange. In the first invocation it just Hence, ๐‘›๐‘š๐‘Ž๐‘ฅโˆ’๐‘…๐‘‡ ๐‘‹3090 < 1.5 G. 32 bit numbers are sufficient to store contains 1 cell. In the second invocation 3, and the ๐‘› column, 2๐‘› + 1. the pyramid values as ED cannot be bigger than the string length. In the previous approach it seems a waste of resources to have most Ignoring the amount of memory of the GPU, our implementation does of the ๐‘š + ๐‘› + 1 workitems working in useless computations. With this not support strings longer than 232 . Longer strings would require a change, the test-case execution time is reduced from 65.85 to 54.24 s, piece-wise strategy and the use a bigger integer types. which represents 190 GCUPS and a speedup factor of almost 1.21ร— with respect to the baseline. 4.1.1. Baseline We could expect that the execution time of the baseline is propor- In a first baseline version, we implement a kernel in which each tional to ๐‘‘ ร— (๐‘š + ๐‘› + 1), where d is the ED between ๐‘ƒ and ๐‘‡ . We workitem computes only one cell from the column. Every kernel in- also could expect that the dynamic workitems improvement execution vocation computes a column from the pyramid. The algorithm has a time is proportional to (๐‘‘ + 1)2 . Using this estimations we compute the dependency between columns, so we must wait until a whole column expected speedup factor from our files, which is 4.2 ร—. But, in reality, 5 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Fig. 6. Conceptual time-line of the interaction between host and OpenCL kernel ฬ‚ โ€™s columns on the motivating example as a depending on the number of final condition checks. Blue areas depict kernel execution Fig. 5. Kernel execution time to compute ๐‘Š in the GPU. Red areas depict the host code manages kernel invocation, while orange function of column number (radius) for 3 different kernel invocation strategies. Baseline areas depict the host code that determines the final condition of the algorithm. Arrows and dynamic strategies suffer a constant overhead caused by the additional time used depict memory transfers in both directions (host-to-device, and device-to-host). (a) doing unnecessary memory transfers to check for the final condition. In addition, the Baseline, where final condition is checked after every kernel invocation. (b) Reducing baseline shows a greater execution time than the dynamic version because it always the number of final condition checks reduces the number of memory transactions and creates more workitems than the other. eliminates the delay between successive kernel executions ๐‘‘๐‘–๐‘˜ . we observe that execution time is not exactly linear with the number of workitems as radius progresses (see cyan line in Fig. 5). What happens is that read operations dominate over kernel execution slowing the overall execution. However, we can see that, in the worst execution time regime, for the same number of active cells the baseline kernel still takes more time to execute. This is because it always launches the worst case distance number of workitems, while the dynamic always launches less. 4.1.3. Minimizing reads To check if a solution is found, the host must check if the furthest reaching radius computed by the kernel on the solution diagonal al- ready reached the final point. This requires to perform an OpenCL read transaction to transfer memory from the device to the host. As depicted in Fig. 6(a), the transaction time can be of the same order, or even longer, than the kernel execution time, becoming a bottleneck. In addition, the memory transactions introduce a delay ๐‘‘๐‘–๐‘˜ that prevents Fig. 7. Effect of reducing the number of read transactions on performance for various the streaming execution of successive kernel invocations. string lengths. Reducing the number of read transactions per kernel invocation below a To reduce the impact of read transactions, we can enqueue several factor 1/100 does not provides any significant additional benefit for all string lengths. kernel invocations before reading the furthest reaching value in the solution diagonal. This strategy increases the number of executed invo- cations by a constant value but this becomes irrelevant when dealing knowledge of the algorithm. In the following subsections we analyze with very long strings. Fig. 6(b) illustrates how the execution time is some strategies to try to increase data locality and, as a result, achieve reduced by grouping kernel invocations before a read operation. better performance. When using this strategy in the example, the execution time is reduced from 55.53 to 37.18 s, which is a 1.49ร— factor with respect 4.2.1. Pre-fetching to shared memory to the execution time of the previous optimization. The performance is As illustrated in Fig. 3, the processing of a pyramid cell requires increased to 277 GCUPS with a workgroup size of 128 workitems. This to access three cells from the previous pyramid column. Hence, each has an impact, since the GPU is active more time without waiting for computed cell (which is stored in global memory) is accessed 3 times the host to interpret the data. The green line in Fig. 5 shows how now the execution time of the kernel basically depends on the number of when computing the following column. elements of each column. Pre-fetching is an optimization strategy [31] that involves copying We analyze how the performance is affected by the reduction of read necessary data from a slower, distant memory to a faster, closer mem- operations with several string lengths (see Fig. 7). After the analysis we ory. When multiple accesses are made to the same data, the overhead decide that 1 read transaction every 100 kernel invocations is a good of the initial copy is offset by the faster subsequent memory accesses on tradeoff to get a good performance on various string lengths. Higher the faster memory. In our case, we try to reduce the memory bandwidth values do not provide additional gains. to global memory by using an strategy based on pre-fetching data from global memory to shared memory. Fig. 8 depicts this strategy. Each 4.2. Strategies to reduce access to global memory workitem works with a row of the pyramid and only fetches data from the cell in the previous column. But for each cell, the algorithm needs to From the programmer perspective, previous designs access the pyra- access 2 additional neighboring cells. The benefit of pre-fetching is that mid ๐‘Šฬ‚ with read and write transactions to global memory. The GPUโ€™s these additional accesses are done on shared memory rather than global cache hierarchy effectively prevents to pay the price of long latencies memory. A corner case happens in the workitems on the boundaries of to global memory, as many accesses result in cache hits. However, the workgroup. They have to do an extra access to global memory since the cache hierarchy acts in a transparent way without an specific one of the neighboring cells is not computed in the workgroup. 6 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Fig. 8. Strategy to fetch elements from global memory to shared memory before computing pyramid cells. Workitems in a workgroup can reuse the data fetched by their sibling workitems after being sure that everyone completed the fetch with a synchronization barrier (illustrated by a blue vertical line). Workitems at the boundaries of the workgroup require an extra access to global memory. Fig. 9. Example of a tiling approach with ๐‘ก = 3 where a workitem computes all cells in a tile. Colored cells depict the cells that must be shared with the workitems in the This approach reduces the memory accesses to global memory next vertical band. by three but requires to synchronize workitems in the workgroup. After fetching their related cell from global memory, all workitems execute a barrier instruction (syncthreads in CUDA) to ensure band, the number of addressed columns is 2๐‘ก, which is higher than the that neighboring workitems in the workgroup already fetched their 2 columns addressed in previous approaches. At a first sight this could data. Workitems at the boundary of the workgroup must access global impose a significantly higher memory demand. However, we notice memory, since one of the neighboring cells is part of another work- that only two columns following a zig-zagged pattern are effectively group. accessed. In Fig. 9, we highlight the ziz-zag data dependencies from Fig. 10(a) shows the that this strategy is counterproductive. The ex- previous tiles with blue and yellow colors. The kernel prefetchs all the pected benefit from reducing accesses to global memory does not occur depending cells from the previous tile band boundary, computes the tile because of two reasons. First, because in the original version, cells that using the local storage, and stores the boundary cells back into global are accessed from multiple workitems are already present in the cache memory. hierarchy of the GPU, so there is no significant difference between the In GPUs, the obvious local storage is shared memory. The shared global memory bandwidth used by both versions. Secondly, because memory is not very big (128 kB per SM in RTX 3090), but is faster the shared memory version introduces a synchronization primitive that than global memory because it is located on-chip. As multiple threads delays the execution of some workitems. will be executed in the SM (the workgroup), a code using the shared memory should divide its use among threads of the workgroup. This 4.2.2. Tiling means that we must allocate 2๐‘ก2 ร— ๐‘ค๐‘”๐‘  bytes for the tiles, where ๐‘ค๐‘”๐‘  The previous shared memory approach suffers from too much syn- is the workgroup size, and then use the workitem local index to index chronization. A way to reduce the synchronization needs is to increase the tile we have to use. the amount of work that a workitem is processing. Tiling is used in In addition to the use of shared memory, we also investigate two most GPU solutions addressing the SW DP algorithm [25]. The goal of additional options for the local storage: global memory and registers. tiling is to improve data locality. The use of global memory can be interpreted as contradictory with our A tiled approach is designed by identifying the data dependencies goal, but we want to analyze if the memory access pattern inside the of the application and extracting the blocks that, after some initial tile has a beneficial effect in the cache hierarchy fetching policies. The data fetching, can be computed without additional accesses to the use of the register file as a memory store can be somewhat challenging. distant data store. In our case we can divide the pyramid in diamond The register file of a GPU is a big memory (compared with classic shaped tiles that, after fetching data from the tile boundary, can be processors register files) used to store registers associated with threads. computed independently. We define tiles by a tile length ๐‘ก parameter, Every thread has its own set of registers, which allow a swap-free which corresponds to the length of the side of the parallelogram. The context switch between threads. In addition, the number of registers number of cells in the tile is 2๐‘ก2 . Tile length can set arbitrarily, while it of a thread is not defined by the architecture (as in most processors might be limited by the required local memory. Increasing tile length like RISC-V having 32) but dynamically defined when kernel code is increases the work done by a workitem and minimizes accesses to compiled. As a result, the register file is divided by the number of global memory. It also results in launching less workitems, less kernel registers required by a thread to know the number of simultaneous on- invocations, and an increase the computation/communication factor, flight threads in a SM. In RTX 3090 every block has 65 536 registers. which is usually beneficial to increase throughput. The use of registers to store the local tile memory is more intricate as When using tiling, compared with the non-tiled approach, the host it requires to describe all tile cells as variables in the source code and must do less invocations and each invocation must use less workitems. implementing additional addressing functions to write and read them. Each kernel invocation computes a vertical tile band instead of the The results of these analyses are shown in Figs. 10(b), 10(c), 10(d). column computed in previous approaches. In the computation of a tile The ๐‘ฅ axis of the figures describe the number of workitems in the 7 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Fig. 10. Performance achieved in NVIDIA GeForce RTX 3090 with different methods that try to reduce access to global memory as a function of the workgroup size. GCUPS are computed as defined in Eq. (2). (a) Pre-fetching in shared memory does not provide a significant with respect to the previous approach due to synchronization overheads. Tiling using one workitem per tile does not provide significant gains as well. Diagrams (b), (c), (d) show multiple tests using different tile lengths (๐‘ก). Performance for tile lengths greater than 3 decreases due to thread divergence and shared memory is the best option for this tiling approach. (e) Having multiple workitems collaborating in a tiled approach provides significant gains. (f) A substantial gain is obtained with the dynamic diamond algorithmic optimization. workgroup. This has an impact on GPU occupancy as higher workgroup the local memory. In these iterations, some workitems will do useless sizes require more memory. Multiple tile lengths are tested, starting work since they do not have anything to compute. Finally, the local from a tile length of 2 increasing tile lengths are drawn with lighter memory is transferred to global memory. The execution of the steps color. It can be observed than generally bigger tile lengths perform is synchronized with the workitems collaborating in the computation worse. done by the workgroup by using a barrier. The use of global memory for local storage does not provide any The โ€˜performance of this design is influenced by the workgroup size, benefit over previous approaches. The achieved best performance is which is determined to the tile size ๐‘ค = 2๐‘ก + 1. This approach provides 256 GCUPS with a tile length of 2 and 128 workitems per workgroup. A a significant gain with respect to previous designs. The maximum profiling analysis reveals that memory accesses are dispersed in mem- achieved performance is 372 GCUPS, when using 381 workitems per ory reducing the chance to coalesce memory operations and increasing workgroup, which corresponds to tile length of 190. the chance to have threads waiting for long memory operations. As a result, the number of non-available threads in the GPU scheduler is 4.3. Dynamic diamond increased. When using registers for the local store, we use conditionals to We left the dynamic diamond approach for the final step, as the implement the equivalent to address decoding. Due to the diamond code to control the elements of the pyramid the must be processed shape of the tile, some diagonals are more frequently accessed than has a higher complexity. When applying this approach, we basically others, so it is convenient that they are checked before less frequent reduce the number of tiles that are computed. This reduction factor ones. A convenient conditional ordering provides a 10 GCUPS gain with is not constant and depends on the distance between both strings. In respect to a naive approach. However the best achieved performance this case, the maximum achieved performance is 516 GCUPS, when with the register-based local store is 267.2 GCUPS with a tile length of using the configuration used in the last iteration (381 workitems per 3 and a workgroup size of 128 workitems. workgroup, which corresponds to tile length of 190). Although not being very significant, the shared memory local store provides an improvement with respect to previous designs. The best achieved performance is 289.6 GCUPS with a tile length of 2 and a 4.4. Summary workgroup size of 224 workitems. After an iterative improvement process, we have achieved a 3.4ร— 4.2.3. Parallel workgroup speed-up factor with respect to the initial WFA design. The main In previous tiling strategies, a whole tile was computed by every elements that make it possible to achieve this speedup factor are: an workitem. Another approach is studied, where several workitems to appropriate kernel invocation strategy that minimizes memory trans- collaborate to compute a tile. The principle of this approach is that each actions between the host and the device, the use of tiling to increase workitem within the workgroup is responsible to fetch the memory and data locality within the kernel, and leveraging the diamond-shaped compute the cells of a single diagonal. Fig. 11 depicts the strategy. In a data dependency to reduce the amount of data to compute. first step, all workitems fetch the two dependent positions from global Fig. 12 illustrates the improvements of each step by reporting the memory and store them in a local store using shared memory. Notice maximum achieved performance (in GCUPS) of every method when that this requires 2๐‘ก + 1 workitems (and not 2๐‘ก โˆ’ 1) because the cells computing the distance between the strings of the illustrating example. in the extremes of the tile access cells from upper and lower diagonals. The presented implementation is open-source, and its code is pub- Then, a loop of 2๐‘ก iterations is required to compute the tile values using licly available at https://github.com/davidcastells/WavefrontVariants. 8 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Fig. 11. Parallel computation of a tile of ๐‘ก = 3. In a first step values are pre-fetched from global memory to shared memory. The computation of tile columns is synchronized with barrier primitives (depicted in blue). Table 2 Execution time of this work on a NVIDIA GeForce RTX 3090 compared with the results from Sadiq and Yousaf [23] (2020) in a cluster running an hybrid MPI OpenMP application. Exp. String A Size String B Size Distance Similarity Reference This work Execution time (s) GCUPS Execution time (s) GCUPS Speed-up factor 1 gbgss201 156,931 gbpln104 79,314 91,590 41% 42.44 0.29 0.12 100.22 353.7ร— 2 gbgss201 156,931 gbhtg11 606,452 450,982 25% 83.44 1.14 0.86 109.92 97.0ร— 3 gbhtg11 606,452 gbgss116 1,517,819 969,770 36% 606.71 1.51 4.81 191.11 126.1ร— 4 gbuna1 308,453 gbinv32 3,424,429 3,116,517 5% 553.04 1.90 6.04 111.83 91.5ร— Table 3 Execution time of this work on a NVIDIA GeForce RTX 3090 compared with the first-phase of MASA-Cudalign on the same platform. Exp. String A Size String B Size Distance Similarity Reference This work Execution GCUPS Execution GCUPS Speed-up time (s) time (s) factor 1 NC_000898.1 162,114 NC_007605.1 171,823 92,254 46% 0.05 557.10 0.13 202.50 0.4ร— 2 NC_003064.2 542,868 NC_000914.1 536,165 271,907 49% 0.55 529.21 0.79 366.20 0.7ร— 3 CP000051.1 1,044,459 AE002160.2 1,072,950 549,784 48% 2.00 560.33 2.73 409.05 0.7ร— 4 BA000035.2 3,147,090 BX927147.1 3,282,708 1,528,200 53% 20.00 516.55 20.81 496.39 1.0ร— 5 AE016879.1 5,227,293 AE017225.1 5,228,663 2,595 99% 36 759.22 0.34 78,734.09 103.7ร— 6 NC_005027.1 7,145,576 NC_003997.3 5,227,293 3,700,051 48% 78.00 478.87 103.38 361.28 0.8ร— 7 NT_033779.4 23,011,544 NT_037436.3 24,543,557 12,565,054 48% 1,200.00 470.65 1,310.61 430.93 0.9ร— 8 BA000046.3 32,799,110 NC_000021.7 46,944,323 15,700,576 66% 2,258.00 681.90 2,159.39 713.03 1.0ร— chr1 NC_000001.10 249,000,000 NC_006468.3 228,000,000 50,882,030 79% 92,162.00 616.00 32,060.21 1,775.17 2.8ร— chr19 NC_000019.9 59,128,983 NC_006486.3 63,644,993 15,090,152 76% 6,557.00 573.93 2,779.20 1,354.08 2.3ร— chr20 NC_000020.10 63,025,520 NC_006487.3 61,729,293 7,609,553 87% 5,933.00 655.74 889.55 4 373.53 6.6ร— chr21 NC_000021.8 48,129,895 NC_006488.2 46,489,110 4,000,655 91% 3,290.00 680.10 167.85 13,330.27 19.8ร— chr22 NC_000022.10 51,304,566 NC_006489.3 49,737,984 6,818,507 86% 3,827.00 666.78 681.13 3,746.37 5.6ร— 23M NT_033779.4 23,011,544 NT_037436.3 24,543,557 12,565,054 48% 1,196.00 472.23 1,309.54 431.28 0.9ร— 10M NC_017186.1 10,236,779 NC_014318.1 10,236,715 374 99.9% 136.00 770.52 0.65 161,216.91 209.2ร— 5M AE016879.1 5,227,293 AE017225.1 5,228,663 2,595 99% 36.00 759.22 0.33 82,823.50 109.0ร— 5. Results In this section, we evaluate our best design and compare it with the state of the art. Since very few works are addressing ED on very long strings we compare our results with the only work fulfilling this requirement and with the first step of the MASA-Cudalign [29] which has a complexity similar to ED. To base our experiments, we use the test files used by [23] and the CUDAlign project publications [26,28] which cover enough differences in terms of string similarity and length. In the case of [23] the files are not exactly the same reported in the original paper as they were the result of concatenating multiple genomic sequences combined in a .seq file, and it has been evolving over time. We have repeated the process with the same sequences files and similar string lengths. The results are similar but not exactly the same reported in the paper. Fig. 12. Maximum performance achieved by the analyzed approaches on the motivating example. The execution platform consist of a Intelยฎ Xeonยฎ W-2155 CPU with 10 cores running at 3.30 GHz, and a NVIDIAโ„ข GPU GeForce RTX 3090, with 82 SMs running at 1.6 GHz. All GPU measurements are using a single GPU. 9 D. Castells-Rufas Parallel Computing 116 (2023) 103019 Table 4 Maximum length and performance of works in the literature computing the ED between two strings. Work Method/Platform Max reported length Max performance (GCUPS) [22] (2021) Banded-BPM/FPGA 300 5076.00 [8] (2022) WFA/GPU 1000 (90% sim.) 10,060.30 (98% sim.) 25,062.60 [15] (2014) BPM/GPU 1000 2300.00 [32] (2018) BPM/CPU ? 1620.00 [17] (2017) ED/GPU 22,500 0.88 [20] (2015) ED/GPU 128,000 20.48 [23] (2020) ED/CPU 3,424,429 1.90 This work WFA/GPU 249,000,000 (41% sim.) 100.22 (79% sim.) 1775.17 (99.9% sim.) 161.216.91 6. Conclusions In this paper, we have presented a GPU implementation of the WFA algorithm adapted to the computation of the ED problem for very long strings. To the best of our knowledge this is the first work considering the problem of computing the ED for strings longer than one hundred million characters. Our baseline GPU implementation is already orders of magnitude faster than previously reported alterna- tives specifically targeting ED. In addition, we present a tiled approach together with the dynamic diamond optimization that contribute to increase the performance by a factor of 3.4ร— with respect to our baseline implementation. For low similarity strings MASA-Cudalign results suggest that its adaptation to the ED problem would provide a competitive perfor- mance. However, for long strings with similarities above 70%, to the best of our knowledge, our results are the best ever reported. In this work we have focused on obtaining the ED value. In future works, we Fig. 13. Our implementation performance (in GCUPS) depends on the similarity of the will focus on obtaining the ED alignment. strings being compared. Declaration of competing interest The Table 2 compares the results of the proposed method in ex- One or more of the authors of this paper have disclosed potential or ecution platform with Sadiq and Yousafโ€™s results reported in the pa- pertinent conflicts of interest, which may include receipt of payment, per [23]. The input strings have a rather low similarity below 50%, but either direct or indirect, institutional support, or association with an Sadiq and Yousafโ€™s top performance is below 2 GCUPS. Our achieved entity in the biomedical field which may be perceived to have poten- speed-up factor (using a single GPU) is always above two orders of tial conflict of interest with this work. For full disclosure statements magnitude with respect to Sadiq and Yousafโ€™s work (which is using refer to https://doi.org/10.1016/j.parco.2023.103019. David Castells- a computing cluster), achieving up to a 353ร— factor for the higher Rufas reports financial support was provided by European Regional similarity file. Development Fund. For the files used in [26,28] we compare with a compilation in our execution platform of the latest version of MASA-Cudalign soft- Data availability ware available at https://github.com/edanssandes/MASA-CUDAlign. The original papers were reporting a performance between 20 and The data used is available from public repositories (data sources 30 GCUPS using a single GPU. But, as shown in Table 3, our tests have been already described in the paper) reveal that the performance in latest GPUs goes up to the region of 600 GCUPS. In this case, our implementation performs slightly worse when Acknowledgments similarity is low, as our performance for similarities lower than 70% falls in the range of 200 to 700 GCUPS. However, for higher similarity This research was supported by the European Union Regional De- strings we easily jump above the TCUPS range reaching a top value of velopment Fund (ERDF) within the framework of the ERDF Operational 161 TCUPS for almost identical files. This dependence of performance Program of Catalonia 2014โ€“2020 with a grant of 50% of the total on similarity is depicted in Fig. 13, which collects all the results from cost eligible under the Designing RISC-V based Accelerators for next the table. String length also has some impact on performance. For the generation computers project (DRAC) [001-P-001723], in part by the same similarity value, longer sequences will generally provide higher Catalan Government under grant 2017-SGR-1624, and in part by the throughput. Spanish Ministry of Science, Innovation and Universities under grant We compare our obtained performance with previous works on the RTI2018-095209-B-C22. literature specifically addressing ED. Table 4 shows the results sorted by the maximum reported length. For short strings (below thousands References of characters) a high performance in the order of TCUPS is generally [1] V.I. Levenshtein, et al., Binary codes capable of correcting deletions, insertions, achieved by processing multiple strings concurrently. In this case, and reversals, Sov. Phys. Dokl. 10 (8) (1966) 707โ€“710. the problem becomes embarrassingly parallel and the performance is [2] J. Wang, Y. Dong, Measurement of text similarity: a survey, Information 11 (9) limited by the communication bandwidth to access input data. (2020) 421, http://dx.doi.org/10.3390/info11090421. 10 D. Castells-Rufas Parallel Computing 116 (2023) 103019 [3] F. Foscarin, F. Jacquemard, R. Fournier-Sโ€™niehotta, A diff procedure for music [19] A. Tomiyama, R. Suda, Automatic parameter optimization for edit distance score files, in: 6th International Conference on Digital Libraries for Musicology, algorithm on GPU, in: International Conference on High Performance Computing 2019, pp. 58โ€“64, http://dx.doi.org/10.1145/3358664.3358671. for Computational Science, Springer, 2012, pp. 420โ€“434, http://dx.doi.org/10. [4] M. ล oลกiฤ‡, M. ล ikiฤ‡, Edlib: a C/C++ library for fast, exact sequence alignment 1007/978-3-642-38718-0_38. using edit distance, Bioinformatics 33 (9) (2017) 1394โ€“1395, http://dx.doi.org/ [20] M. Kruliลก, D. Bednรกrek, M. Brabec, Improving parallel processing of matrix-based 10.1093/bioinformatics/btw753. similarity measures on modern gpus, in: International Conference on Similarity [5] S. Arora, B. Barak, Computational Complexity: A Modern Approach, Cambridge Search and Applications, Springer, 2015, pp. 283โ€“294, http://dx.doi.org/10. University Press, 2009. 1007/978-3-319-25087-8_27. [6] R.A. Wagner, M.J. Fischer, The string-to-string correction problem, J. ACM 21 [21] D. Bednรกrek, M. Brabec, M. Kruliลก, Improving matrix-based dynamic program- (1) (1974) 168โ€“173, http://dx.doi.org/10.1145/321796.321811. ming on massively parallel accelerators, Inf. Syst. 64 (2017) 175โ€“193, http: [7] S. Marco-Sola, J.C. Moure, M. Moreto, A. Espinosa, Fast gap-affine pairwise //dx.doi.org/10.1016/j.is.2016.06.001. alignment using the wavefront algorithm, Bioinformatics 37 (4) (2021) 456โ€“463, [22] D. Castells-Rufas, S. Marco-Sola, Q. Aguado-Puig, A. Espinosa-Morales, J.C. http://dx.doi.org/10.1093/bioinformatics/btaa777. Moure, L. Alvarez, M. Moretรณ, OpenCL-based FPGA accelerator for semi-global [8] Q. Aguado-Puig, S. Marco-Sola, J.C. Moure, D. Castells-Rufas, L. Alvarez, A. approximate string matching using diagonal bit-vectors, in: 2021 31st Interna- Espinosa, M. Moreto, Accelerating edit-distance sequence alignment on GPU tional Conference on Field-Programmable Logic and Applications (FPL), IEEE, using the wavefront algorithm, IEEE Access 10 (2022) 63782โ€“63796, http: 2021, pp. 174โ€“178, http://dx.doi.org/10.1109/FPL53798.2021.00036. //dx.doi.org/10.1109/ACCESS.2022.3182714. [23] M.U. Sadiq, M.M. Yousaf, Distributed algorithm for parallel edit distance com- [9] K.E. Iverson, A programming language, in: Proceedings of the May 1-3, 1962, putation, Comput. Inform. 39 (4) (2020) 757โ€“779, http://dx.doi.org/10.31577/ Spring Joint Computer Conference, 1962, pp. 345โ€“351. cai_2020_4_757. [10] D.S. Hirschberg, A linear space algorithm for computing maximal common [24] T.F. Smith, M.S. Waterman, et al., Identification of common molecular subse- subsequences, Commun. ACM 18 (6) (1975) 341โ€“343, http://dx.doi.org/10. quences, J. Mol. Biol. 147 (1) (1981) 195โ€“197, http://dx.doi.org/10.1016/0022- 1145/360825.360861. 2836(81)90087-5. [11] G. Myers, A fast bit-vector algorithm for approximate string matching based on [25] Z. Xia, Y. Cui, A. Zhang, T. Tang, L. Peng, C. Huang, C. Yang, X. Liao, A review dynamic programming, J. ACM 46 (3) (1999) 395โ€“415, http://dx.doi.org/10. of parallel implementations for the Smithโ€“Waterman algorithm, Interdiscip. Sci.: 1145/316542.316550. Comput. Life Sci. (2021) 1โ€“14, http://dx.doi.org/10.1007/s12539-021-00473-0. [12] E. Ukkonen, Algorithms for approximate string matching, Inf. Control 64 (1โ€“3) [26] E.F.O. Sandes, A.C.M. de Melo, CUDAlign: using GPU to accelerate the com- (1985) 100โ€“118, http://dx.doi.org/10.1016/S0019-9958(85)80046-2. parison of megabase genomic sequences, in: Proceedings of the 15th ACM [13] H. Hyyrรถ, A bit-vector algorithm for computing levenshtein and damerau edit SIGPLAN Symposium on Principles and Practice of Parallel Programming, 2010, distances, Nordic J. Comput. 10 (1) (2003) 29โ€“39. pp. 137โ€“146, http://dx.doi.org/10.1145/1693453.1693473. [14] J. Zhang, H. Lan, Y. Chan, Y. Shang, B. Schmidt, W. Liu, BGSA: a bit-parallel [27] E.F.d.O. Sandes, A.C.M. de Melo, Retrieving smith-waterman alignments with global sequence alignment toolkit for multi-core and many-core architec- optimizations for megabase biological sequences using GPU, IEEE Trans. Parallel tures, Bioinformatics 35 (13) (2019) 2306โ€“2308, http://dx.doi.org/10.1093/ Distrib. Syst. 24 (5) (2012) 1009โ€“1021, http://dx.doi.org/10.1109/TPDS.2012. bioinformatics/bty930. 194. [15] A. Chacรณn, S. Marco-Sola, A. Espinosa, P. Ribeca, J.C. Moure, Thread- [28] F.d.O. Edans, G. Miranda, A.C. de Melo, X. Martorell, E. Ayguadรฉ, Cudalign cooperative, bit-parallel computation of levenshtein distance on GPU, in: 3.0: Parallel biological sequence comparison in large gpu clusters, in: 2014 Proceedings of the 28th ACM International Conference on Supercomputing, 2014, 14th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing, pp. 103โ€“112, http://dx.doi.org/10.1145/2597652.2597677. IEEE, 2014, pp. 160โ€“169, http://dx.doi.org/10.1109/CCGrid.2014.18. [16] K. Balhaf, M.A. Shehab, T. Walaโ€™a, M. Al-Ayyoub, M. Al-Saleh, Y. Jararweh, [29] E.F. de Oliveira Sandes, G. Miranda, X. Martorell, E. Ayguade, G. Teodoro, A.C.M. Using gpus to speed-up levenshtein edit distance computation, in: 2016 7th Melo, CUDAlign 4.0: Incremental speculative traceback for exact chromosome- International Conference on Information and Communication Systems (ICICS), wide alignment in GPU clusters, IEEE Trans. Parallel Distrib. Syst. 27 (10) (2016) IEEE, 2016, pp. 80โ€“84, http://dx.doi.org/10.1109/IACS.2016.7476090. 2838โ€“2850, http://dx.doi.org/10.1109/TPDS.2016.2515597. [17] K. Balhaf, M.A. Alsmirat, M. Al-Ayyoub, Y. Jararweh, M.A. Shehab, Accelerating [30] T. Allen, R. Ge, Characterizing power and performance of gpu memory access, levenshtein and damerau edit distance algorithms using GPU with unified mem- in: 2016 4th International Workshop on Energy Efficient Supercomputing (E2SC), ory, in: 2017 8th International Conference on Information and Communication IEEE, 2016, pp. 46โ€“53, http://dx.doi.org/10.1109/E2SC.2016.012. Systems (ICICS), IEEE, 2017, pp. 7โ€“11, http://dx.doi.org/10.1109/IACS.2017. [31] P. Hijma, S. Heldens, A. Sclocco, B. van Werkhoven, H.E. Bal, Optimization 7921937. techniques for GPU programming, ACM Comput. Surv. (2022) http://dx.doi.org/ [18] Y. Li, L. Schwiebert, Memory-optimized wavefront parallelism on GPUs, Int. J. 10.1145/3570638. Parallel Program. 48 (6) (2020) 1008โ€“1031, http://dx.doi.org/10.1007/s10766- [32] Y. Chan, K. Xu, H. Lan, B. Schmidt, S. Peng, W. Liu, Myphi: efficient levenshtein 020-00658-y. distance computation on xeon phi based architectures, Curr. Bioinform. 13 (5) (2018) 479โ€“486, http://dx.doi.org/10.2174/1574893612666171122150933. 11