GPU-accelerated CBC search

Acceleration of low-latency gravitational wave searches using Maxwell-microarchitecture GPUs


Low-latency detections of gravitational waves (GWs) are crucial to enable prompt follow-up observations to astrophysical transients by conventional telescopes. We have developed a low-latency pipeline using a technique called Summed Parallel Infinite Impulse Response (SPIIR) filtering, realized by a Graphic Processing Unit (GPU). In this paper, we exploit the new Maxwell memory access architecture in NVIDIA GPUs, namely the read-only data cache, warp-shuffle, and cross-warp atomic techniques. We report a 3-fold speed-up over our previous implementation of this filtering technique. To tackle SPIIR with relatively few filters, we develop a new GPU thread configuration with a nearly 10-fold speedup. In addition, we implement a multi-rate scheme of SPIIR filtering using Maxwell GPUs. We achieve more than 100-fold speed-up over a single core CPU for the multi-rate filtering scheme. This results in an overall of 21-fold CPU usage reduction for the entire SPIIR pipeline.

04.80.Nn, 95.75.-z, 97.80.-d, 97.60.Gb

1 Introduction

We are entering an exciting time of gravitational wave (GW) astronomy. The first GW signal is detected in September 2015 by the Laser Interferometric Gravitational-Wave Observatory (LIGO) [1]. This opens up a new window to multi messenger astronomy with unprecedented power of discovery, in probes of some of the most enigmatic transients in the sky for their emissions in both the electro-magnetic and gravitational wave spectrum, e.g., short or long gamma-ray bursts produced in binary coalescence and core-collapse supernovae [2, 3]. In particular, low-latency detection and localization of GW sources are gaining priority in order to enable prompt electromagnetic (EM) follow-up observations of GW sources. Capture of transient EM events triggered by low-latency GW events will bring breakthroughs in understanding the mechanism of sources [4].

The optimal method to search for compact binary coalescence (CBC) GW sources is by correlating the expected waveform templates with data, known as the matched filtering method [5, 6]. An efficient way to compute the correlation is by the fast Fourier transform (FFT) technique, especially so when the array size is large. Low latency searches, however, require real-time processing of relatively short data-segments, limiting any gain obtained in FFT-based correlations. Existing efforts of low-latency GW detection from binary coalescence sources by LIGO-Virgo community employ sophisticated strategies to realize low-latency with manageable computation, notably a Multi-Band Template Analysis (MBTA) [7] method utilizing the two-band frequency-domain Finite-Impulse-Response (FIR) filters to represent a given template for filtering , the Low-Latency Online Inspiral Detector (LLOID, a.k.a GstLAL for the software name) [8] method utilizing multi-band FIR filters and the Summed Parallel Infinite Impulse Response (SPIIR) [9, 10, 11] method utilizing infinite impulse response (IIR) technique in the time domain which is expected to have zero latency. These three techniques evolve into automated software programs, known as pipelines, which additionally allow input from multiple detectors and employ template-based statistical tests to veto transient noises and produce GW alerts to the community. As of this writing, these three pipelines have achieve a medium latency of less than one minute. This paper is focused on the SPIIR method. Compared to the FFT technique, it is expected to be more efficient at latencies lower than tens of seconds for advanced detectors [9].

Our previous work of GPU accelerated SPIIR filtering method uses NVIDIA’s Fermi GPUs with a speedup a factor in the order of 50-fold over a single core Intel i7 CPU [12]. However, that GPU optimization targeted the SPIIR filtering with number of SPIIR filters in the range of 128 to 256. While in the multi-rate filtering scheme, the total filters are split to be applied to data at different sampling rates, that the number of filters can be as small as a few in some sampling rates. In this paper, we extend the GPU optimization of SPIIR filtering for various number of filters and explore the features of the recent Maxwell GPUs. Our optimization here achieves 3-fold improvement over the previous targeted range and achieves up to 10-fold speedup beyond the range, compared to the previous work. Compared to the SPIIR filtering on a single-core CPU, our GPU acceleration is now 60 to 125 fold faster depending on the number of filters to be applied.

In addition to the SPIIR filtering, we have extended the GPU acceleration to other components of the SPIIR pipeline. Another bottleneck of the SPIIR pipeline in terms of computational efficiency is the sampling-rate alterations of multi-rate filtering. In multi-rate filtering, a set of filters is applied to data at different sampling rates, the results of which are combined to the initial rate. Our acceleration of the multi-rate filtering scheme realizes a -fold improvement in CPU resource reduction and hence a fold reduction in CPU resource for the entire SPIIR pipeline. The filtering process in this paper uses single-precision floating-point number format, the result of which has ignorable difference with the results using double-precision format.

The structure of this paper is as follows: In Sec. 2 we review the mathematical representation of the SPIIR method and present the multi-rate filtering scheme for our GPU implementation. In Sec. 3 we present the technical details of our new GPU optimization on SPIIR filtering exploiting the newest features from Maxwell GPUs and GPU optimization technique on sampling-rate alterations. In Sec. 4 we present the performance results of our GPU accelerated SPIIR filtering and the multi-rate filtering scheme. In Sec. 5, we present the performance of our GPU accelerated low-latency SPIIR detection pipeline. The conclusion and future work will be given in Sec.6.

2 Multi-rate SPIIR filtering

2.1 SPIIR method

Matched filtering is known to be the optimal detection method in deep searches for signals in stationary and Gaussian noise. Here, we consider a CBC waveform template in the time domain . The optimal detection output, known as the signal-to-noise ratio (SNR), is given by:


i.e. the cross-correlation between whitened waveform template and the whitened detector strain , which is given by:


where is the Fourier transform of detector data s() and is the one-sided noise spectral density of a detector defined through the expectation :


By sampling at discrete instances at a sampling rate , Eq. 1 can be rewritten in a discrete form,


CBC searches usually adopt the matched filter , that takes Eq. 1 into a convolution integral. Quite different from the common FIR method, SPIIR method utilizes IIR filters to reconstruct the matched filter. A first-order IIR filter bears a simple form shown in Eq. 5 and Eq. 6:


where is the filter output at time step (), is the filter input, and and are complex coefficients. A solution to this first-order linear inhomogeneous difference equation is:


For a target as complex as a CBC signal, a group of first-order IIR filters are developed with each filter representing a small segment of the matched filter [9, 11]. The number of IIR filters that are needed to reconstruct a highly accurate matched filter varies from a dozen to several hundreds, depending on the complexity of the waveform and the the limit of the detection band. The filter construction procedure can be found in [9, 11]. Here, we express the output of a SPIIR filter as:


where is the output for the th filter and is the time-delay for this filter. The discrete form of SNR output from a group of SPIIR filters is given by:


2.2 Multi-rate implementation of SPIIR filtering

Figure 1: Schematic diagram of multi-rate implementation of SPIIR filtering scheme. The input and output sampling rates are Hz. represent the number of filters to be applied to the corresponding rates.

Our implementation of the multi-rate filtering scheme for GPU acceleration is shown in Fig. 1. For implementation conveniency, the data is downsampled by a factor of 2 in succession. Each sample rate stream is filtered using the corresponding SPIIR filters. The filtering output of the lowerest rate will be upsampled and added to the filtering output of next highest rate. This process will repeat several times until reaching the initial sampling rate.

Which sampling rate should a given SPIIR filter work on is essentially determined by the coefficient of the filter. The coefficient determines the upper bound of the frequency band of the filter, and thus the Nyquist rate [13] for the filter to function. We round the Nyquist rate of the filter to the nearest biggest available rate for the filter to work on.


To avoid the known problems of spectral leakage caused by the squared window, we adopt the popular Kaiser-window-tapered low-pass filter implemented in the open-source gstreamer library 4 to perform the interpolation for resampling. The resampling formula is shown in Eq. 9 where is the data, and represents the Kaiser-window-tapered low-pass filter. and are discrete sampling points. is the original sampling rate, is the resampled rate. is assumed to be bandlimited to . is the length of the filter. We call the Kaiser-window-tapered low-pass filter ‘Kaiser filter’ in the rest of the paper.

There is a single parameter that controls the quality, measured by the stop-band attenuation, of Kaiser filter. The gstreamer library provides 11 Kaiser filters with stop-band attenuation up to 100 units of decibels (dB). To avoid the band aliasing of downsampling, we choose the Kaiser filter with stop-band attenuation of dB. We choose the Kaiser filter with stop-band attenuation of dB for upsampling that we found work most efficiently while maintaing signal recovery quality in practice.

The expected computational efficiency of our multi-rate scheme can be estimated as follows. We denote the number of total SPIIR filters of a given template as ; the full-rate we are considering as ; and the number of search templates as . According to Eq. 7 and Eq. 8, filtering on one data point requires 12 floating operations. Thus the total floating point operations per second (flops) for SPIIR filtering at full-rate is . The dB downsampling Kaiser filter has 384 steps in gstreamer and the dB upsampling Kaiser filter has 32 steps. The resampling and summation of filtering result at different sample rates are negligible. If half of the filters can be applied to sub-rate , the cost will reduce by . A factor of a few savings on the computation cost is expected if more filters can be applied to even lower rates..

3 Optimization of Multi-rate SPIIR filtering on Maxwell GPUs

There are several interfaces to use GPU for general purpose problem, including NVIDIA’s Compute Unified Device Architecture (CUDA) [14] language, Khronos Group’s Open Computing Language (OpenCL) [15], Microsoft’s C++ Accelerated Massive Parallelism (C++ AMP) [16]. We use the popular CUDA language for our GPU acceleration.

There are several elements involved in any GPU acceleration, here using CUDA, namely the mapping of algorithm functions to CUDA threads and blocks, and the mapping of data to GPU hierarchical memory architecture. We focus on these optimization strategies accordingly (Sec. 3.1.1, Sec. 3.1.2, Sec. 3.2). In this process, we set out to exploit the advanced GPU memory exchange mechanism of Maxwell GPUs, namely the warp-shuffle and atomic operation techniques (sec 3.1.3).

We provide a general explanation of the relation of a GPU hardware and the CUDA semantics for reference here. A GPU chip consists of several Streaming Multiprocessors (SMs). The Maxwell GPU features improved SM architecture renamed as SMM. One SMM has many processing cores upon which one is capable of create, managing, scheduling, and executing CUDA threads. CUDA threads are executed in groups of 32, which are named ”warps”. A group of CUDA warps aggregate to a CUDA block. While one block is limited to one SMM, one SMM can have several blocks running concurrently. A GPU has a hierarchy of memory spanning a range of access speeds and storage sizes. The choice of memory to use may greatly affect the overall GPU performance.

3.1 Optimization of SPIIR filtering in single-rate

Optimization for varying number of filters

Our previous GPU optimization [12] considers SPIIR filtering with the number of filters a few hundred. It is not highly-optimized toward filtering with small number of filters, which is likely the case in the multi-rate filtering scheme. Here, we develop a CUDA kernel for templates with filters, which we call warp-based kernel, where the filters in one or multiple templates are mapped to one CUDA warp. This allows us to not only avoid branch execution paths in one warp but also significantly reduce the number of idle threads. We keep the block-based CUDA kernel, presented in the previous work, for templates with filters where each template is mapped to one CUDA block with multiple CUDA warps. Fig. 2 shows the actual mapping between IIR filters for one template and GPU threads for one block regarding the two kinds of configurations.

Figure 2: Schematic of the SPIIR template-filters hierarchy mapped onto the CUDA block-warp-threads hierarchy. ‘Mul’ means ‘Multiple’ in this figure. is the number of filters of any given SPIIR template.

The details of mapping of number of filters to CUDA threads is given here. We exploit the fact that CUDA threads are executed in groups of 32 (warp size), launching as much as possible our CUDA kernels with a number of threads = multiple of warp size. This is considered to be optimal as it helps avoid idled threads. For a template with filters, we assign threads for this template, where is rounded to the next multiple of 32 after . For a template with filters, we assigned threads where is rounded to the next power of 2 after . Multiple templates may be executed in a warp if they are able to fit into the warp. For instance, a template with 513 filters will be assigned to 544 threads, while a template with 5 filters will be assigned to 8 CUDA threads. Four small templates (with 5 filters) will be executed within a single warp. With this assignment, it can be guaranteed that the number of idle threads will not be more than 31 in the worst case.

We use the Maxwell architecture-based GeForceGTX 980 GPU as our test machine, it has maximumly 2048 active threads, To maximize the active number of threads in each SM, the minimum number of threads for one block (the number of blocks per SM is 32) must be larger than . For our warp-based CUDA kernel, we chose to use 256 threads for each block as recommended in CUDA Programming Guide [14]. For our block-based CUDA kernel, only one template is mapped to a single block.

Read-only data cache

Data access is a critical aspect which can greatly affect the GPU performance. For unoptimized GPU programs, it may take longer to perform memory access than to do the core calculation. We analyzed the data access pattern of SPIIR filtering and applied an efficient data mapping method to improve the GPU performance. Three types of memory are investigated in our implementation.

Register is the fastest accessible memory. It is local to individual threads and cannot be accessed by other threads. At the block level, there are shared memories accessible by all the threads of the same block (but not necessarily by the same warp). For the GTX 980 GPUs, one Maxwell streaming multiprocessor features 96 KB of dedicated shared memory. Although shared memory features a broader memory access, it is much slower than the register, and usually requires synchronization within the block. Global memory is located off the chip and its speed is the slowest in the GPU memory hierarchy. A new feature is introduced with the Kepler GK110 architecture (a predecessor of the Maxwell architecture)—the read-only data cache. It can be used to cache read-only global memory load which reduces the global memory access time. The GTX 980 has 4GB of global memory shared by all threads. To achieve high global memory throughput, we applied a coalesced memory access technique which can combine multiple memory accesses into as few as one cache transaction.

Figure 3: A schematic on how the data are mapped onto the GPU memory hierarchy to achieve low-latency data access for SPIIR filtering. The left part shows the color-coded GPU memory hierarchy. From top to bottom, the data access speed decreases but the data storage capacity increases. The right diagram illustrates the memory types created for variables for SPIIR filtering. is the input data to be filtered, are filter coefficients, is intermediate output from each filter, represents the iterative process, and SNR is the filtering result.

Fig. 3 is a schematic on how we organized and mapped the data of one SPIIR template to the GPU memory hierarchy in order to achieve low-latency data access. The output of the SPIIR method is given by Eq. 7 and Eq. 8. Due to the iterative nature of the IIR filter, the left-hand-side of Eq. 7, , will be reused by the next iteration. As there are a huge number of iterations, we therefore chose to store into a register to utilize the fastest memory access in GPU. The input parameter cannot be reused by the same filter, but other filters of the same template may need to access it. The access pattern of has a good temporal and spatial locality so they are stored in the CUDA read-only data cache. Other input parameters, such as , cannot be reused and so they are stored in the slower global memory, but we utilized the efficient coalesced memory access feature. Finally, the final SNR outputs are stored into global memory.

Inner-warp and cross-warp summation optimization

Obtaining the final SNR in Eq. 8 requires summation over all results of filters. This suggests a synchronization of outputting individual results when performing parallel computing. Previously we used the implicit synchronization feature of the CUDA warp to significantly reduce the cost of parallel threads synchronization and a multiple-thread parallel sum reduction method to reduce the cost of summation steps [12].

Here the summation of all results within a warp is further improved by the warp-shuffle technique introduced from Kepler michroarchitecture GPUs. The warp-shuffle technique allows threads to read registers from other threads in the same warp. This is a great improvement over the previous higher-latency shared memory exchange within a warp.

For our block-based CUDA kernel, where the template size is larger than 32, the output cannot be calculated without any cross-warp communication. We consider three different cross-warp summation methods to calculate the final SNR from partial SNRs. The first one is the Direct Atomic summation (DA) method. It uses atomic operations of summation in global memory for all the partial SNRs. This is the simplest and straightforward way to calculate SNRs from multiple partial SNRs. Atomic operations were considered expensive and ought to be avoided as much as possible before (without affecting the correctness of results). From Maxwell GPUs, the atomic operations have been improved significantly.

The second method is Shared memory Warp-shuffle (SW) method. It collects all the partial SNRs of N times iterations into the shared memory of one warp (the batched computation model proposed in [12]) and performs the warp-shuffle operation for the final SNR. Therefore, we need additional shared memory operations to calculate the final SNR based on the partial SNRs, one explicit synchronization operation to synchronize all the threads.

The last method is Shared Memory Atomic Summation (SA) method. It collects partial SNRs into shared memory of one warp and performs atomic operation to compute the final SNR. This method also involves the same shared memory loading overhead and synchronization overhead as the SW method. Later we will show that the simplest atomic operations gives the best performance compared with the SW and SA methods.

3.2 Optimization of resampling and summation

Figure 4: Schematic of the downsampling data hierarchy mapped onto the CUDA block-threads hierarchy. ‘Mul’ means ‘Multiple’, and ’1s’ means ’one second’ in this figure.
Figure 5: Schematic of the downsampling GPU kernel function and how data mapped onto the CUDA GPU memory hierarchy. Each color represent a type of GPU memory. Green: shared memory; red: register memory; blue: global memory.
Figure 6: Schematic of the upsampling input data hierarchy mapped onto the CUDA block-threads hierarchy. ‘Mul’ means ‘Multiple’, and ’1s’ means ’one second’ in this figure.
Figure 7: Schematic of the upsampling and combination GPU kernel function and how data mapped onto the CUDA GPU memory hierarchy. The left part shows the GPU memory hierarchy with explanation of the color code in Fig. 5.

We optimize the usage of CUDA threads and various types of memory as in Fig. 5 and Fig. 5 for the downsampling process. As shown in Fig. 5, we map the production of one downsampled data point to one CUDA thread. The number of downsampled data points can vary over a power of 2 series ranging from 32 to 2048. For downsampling with output data points more than the number of GPU cores, the GPU cores will be fully occupied. It is inevitable for some sub-rate downsamplings, there will be idle GPU cores. As downsampling is the least computational process, that it is only performed once for all templates, in our multi-rate scheme, we do not further optimize the idle cores. The number of blocks allocated is designed as where is the number of downsampled points. For the three types of GPU memories, we map the input and downsampling Kaiser filter to the shared memory as they will be reused a number of times. Kaiser filtering is calculated iteratively and the intermediate filtering result is stored in the register memory for fastest data access. Finally, The global memory will store the downsampling output.

Similarly, Fig. 7 and Fig. 7 show our CUDA design for the combined function of upsampling and upstream summation. As a one-second SNR series has a real and an imaginary component, we map each component of the SNR series to one CUDA block. Each thread is mapped to upsampled SNR points where is the number total SNR points in a second. Different from the downsampling that only is perforce once, the upsampling and upstream summation needs to be performed for each template. The number of blocks will double the number of templates. The number of threads will be mostly likely more than the number of GPU cores, making a full occupancy of GPU cores. For the memory mapping, the upsampling Kaiser filter is mapped to the shared memory and the intermediate filtering output is mapped to the register memory, the same as the downsampling memory mapping. The input SNR series can not fit in the shared memory and they are stored in GPU global memory.

4 Results of multi-rate SPIIR filtering

In this section, we first show the GPU performance of SPIIR filtering and the improvement from each optimization step (Sec. 4.1). In the second part (Sec. 4.2)., we show the GPU performance of several scenarios of multi-rate SPIIR filtering. This promises a lucrative performance improvement by GPU acceleration for our SPIIR pipeline.

All GPU implementations are tested on a GeForce GTX 980 (Maxwell microarchitecture) GPU equipped desktop computer where Tab. 1 shows its configuration. As the CPU counterparts are implemented in single CPU thread fashion, we use the elapsed time as our performance measurement to reflect the usage of CPU resource. We run each experiment ten times and use the average time as the timing result. We set the number of templates to 4096 for performance testing purpose in this section. Note that the number of templates for a search can range from a few hundreds to hundreds of thousands.

Hardware CPU Intel Core i7-3770 3.40 GHz
Host Memory 8 GB DDR3
GPU Memory 4 GB DDR5
Software Operating System Fedora 20 64-bit
Host Compilation gcc 4.8.3 -O2
CUDA Compilation nvcc 6.5
Table 1: Testbed configuration

4.1 Performance study of SPIIR filtering in single-rate

In this section , We first show the overall improved performance of our new GPU acceleration (noted as New Kernel) over the previous GPU acceleration [12] (noted as Pre-Kernel). We then show the breakdown of the improvement by exploiting new features of Maxwell GPU.

Varying number of filters

New Kernels and Pre-Kernels are tested against the the filtering using a single-core CPU. Tab. 2 shows the speedup ratio of New Kernels in comparison to Pre-Kernels with different number of filters for filtering. For number of filters where both kernels use the same thread configuration, the improvements by the New Kernels over Pre-Kernels increase as the number of filters increase, as shown by Tab. 2. This is mainly due to that the New Kernels improve the data access and exchange speed, the effects of which are manifested with more number of filters used. For the previous targeted optimization range (i.e. number of filters 128 to 256), the New Kernel have about 3-fold speed improvement over the Pre-Kernel. For the size of template (i.e. the number of filters) is relatively low (template size 32), our new thread configuration — warp-based kernel, improve the speed performance as the number of filters decrease, to a factor nearly 10-fold.

Template Size 4 8 16 32 64 128 256 512
Kernel 63.10 70.62 61.22 55.98 109.76 124.41 125.58 123.97
6.13 10.79 19.61 37.50 42.64 48.72 39.51 29.78
Table 2: Speedup ratios of different GPU kernels compared with the CPU counterpart.

Performance Improvement by Employing Warp-shuffle and Read-Only Data Cache

Fig. 9 shows the performance of two different implementations with and without warp-shuffle. The first implementation uses warp-shuffle to access data within a warp to calculate the summation. The second implementation uses shared memory to access data and calculate summation. It shows that warp-shuffle summation performs significantly faster than shared memory summation.

Fig. 9 shows the effect of using the read-only data cache. The Maxwell GPU loads its global memory data to L2 cache rather than to L1 cache. In Fig. 9, GPU L2 cache illustrates the kernel performance when storing parameter of Eq. 7 in the L2 cache, which is the default cache for global memory. Because variables such as of Eq. 7 has a favorable temporal and spatial locality, using the read-only data cache can be very efficient as shown in the figure.

Figure 8: Performance comparison in performing parallel summation reduction between shared memory and warp-shuffle using 4096 templates

Figure 9: Memory access performance comparison between caching parameter in L2 cache and read-only data cache using 4096 templates

Performance Comparison Among Three Cross-warp Summation Methods

Fig. 10 illustrates the performance of the block-based SPIIR kernels using three different cross-warp summation methods: DA (directly using atomic operation in global memory),SW (using warp-shuffle and shared memory) and SA (using atomic operation in shared memory) methods. To ensure objective comparison, we attempt to optimize each implementation as much as possible. GPU bank conflicts are perfectly avoided in the SW and SA methods where shared memory usage is involved. As one would expect that shared memory access is much faster than global memory access, the SW and SA methods could be faster than DA method. Surprisingly, the DA method surpasses the other two methods as shown in Fig. 10.

Figure 10: Execution time of different cross warp summation methods with 4096 templates in different execution configurations.

We design two additional experiments to figure out the reason for the inferior performance of the shared memory methods (SW and SA), whether it is because of the cost of shared memory accesses or the explicit synchronization. To reduce the total synchronizaiton cost, both SW and SA methods take advantage of shared memory to execute many iterations before one explicit synchronization operation. Fig. 12 clearly illustrates that the cost of synchronization can be significantly reduced or amortized into many iterations. The larger the iteration number, the better the performance. However, when the iteration number reaches a certain point, the performance benefit becomes almost unchanged as the iteration number continue growing. Fig.12 shows that the reason is that the synchronization overhead has been completely hidden by calculation in such circumstances. In Fig.12 we disabled the synchronization operation in these shared memory GPU kernels to observe the performance influence of explicit synchronization operations in sufficiently large iteration number. Though the computation result may not be correct, the comparison itself does show the impact of explicit synchronization operations. As can be seen, the performance between synchronization and no synchronization GPU kernels are almost unidentifiable, illustrating that we can completely remove the overhead of explicit synchronization by improving the iteration number. Therefore, the additional shared memory operations introduced in SW and SA methods lead to their lower performance compared with DA method.

Figure 11: The optimization effect when we change the iteration number. Template size of 256 is used (with 4096 templates).

Figure 12: The overhead of explicit synchronization can be ignored when we use large iteration number. Iteration number 256 is used here (with 4096 templates)

The impact of atomic operations is also tested for the DA method. Atomic operations are typically considered costly and should be avoided whenever possible, however Fig.13 shows the cost of atomic operation can be almostly ignored because the execution time of two kernels with and without atomic operations is so close. All the experiments explain why DA method is better than SW and SA methods.

Figure 13: Execution time of atomic and non-atomic kernels with 4096 templates in different execution configurations.

4.2 Performance of SPIIR filtering in multi-rate

The results shown in the last section is SPIIR filtering in single-rate. Here we show the performance of the GPU-accelerated multi-rate implementation of SPIIR filtering, which includes GPU accelerated rate alterations. We test several multi-rate scenarios denoted as ”number of rates” in Tab. 3. “Number of rates” as 1 in the table denote the single rate at full rate (4096 Hz). The number of filters at full rate for the 4096 templates is set as 1024. The sub-rates are formed according to the design shown in Sec. 2.2. The initial filters are equally divided to each sub-rate for testing purpose. Tab. 3 shows elapsed times for the CPU and the GPU-accelerated multi-rate filtering in different scenarios. The efficiency of our multi-rate design shown here is consistent with our estimation in Sec. 2.2 that the 8-rate scenario reduces the measured time to around one quarter. Our GPU acceleration of multi-rate scheme is dominated by the GPU acceleration of SPIIR filtering in single-rate. Therefore it is very efficient that we are able to obtain more than 100-fold speedup for any scenario.

Table 3: Elapsed times of CPU and GPU-accelerated multi-rate filtering scheme in different scenarios.

5 Results of the GPU-accelerated low-latency SPIIR pipeline

Figure 14: Schematic flowchart of the low-latency SPIIR pipeline gstlal_iir_inspiral with data input from two detectors. Solid blocks depict main components of this pipeline. Dashed blocks depict external processes.

The low-latency SPIIR detection pipeline we are considering here is publicly available through distribution of gstlal software library 5. gstlal provides a variety of components from LIGO Algorithm Library (LAL) 6 for LIGO data processing and uses the gstreamer framework to control streaming data. We use the existing gstlal components for reading calibrated data, data whitening, performing coincident post-processing for our pipeline. We use our own multi-rate SPIIR filtering for template filtering. The pipeline has the code name gstlal_iir_inspiral in gstlal. A schematic flowchart of the pipeline is shown by Fig. 14.

The gstreamer framework inherently employs multi-threading technique to take advantage of the multiple cores of a CPU. Therefore, the CPU implementation of the SPIIR pipeline is by default multi-threaded. We propose a different criterion, the CPU time used by all CPU cores, rather than the elapsed time for the single-threaded CPU implementation, to measure the performance of the GPU accelerated pipeline versus the original CPU pipeline.

Multi-rate SPIIR filtering Other
SPIIR filtering Resampling Upstream summation
93% 3.1% 0.9% 3%
Table 4: Computation profiling of the CPU low-latency SPIIR pipeline gstlal_iir_inspiral.

We present the used CPU time of the main components of the CPU pipeline in percentages, shown by Tab. 4. This computation profiling is performed using the platform given by Tab. 1. The data we process is a short segment of recolored LIGO 5th Science run data, which represents a set of clean data and produces a reasonable number of single events for coincidence analysis. The SPIIR filtering dominates the computation, taking over of the CPU time while the resampling and summation come next at .

We now can estimate the expected CPU resource reduction of the pipeline, if part of the whole the multi-rate filtering component is executed on GPU instead of CPU. The expected resource reduction is determined by the resource reduction brought by accelerated module and the computational cost of this module in the pipeline. It is given by:


If we only apply GPU acceleration only on the module of SPIIR filtering in our pipeline and we choose a somewhat -fold from Tab. 2 for this component. The total resource reduction ratio of our pipeline will be not more than -fold. If on top of that, we apply GPU acceleration on resampling and upstream summation components, and we choose a somewhat -fold from Tab. 3 for this multi-rate SPIIR filtering. The total resource reduction of our pipeline will be about -fold.

We measure the CPU resource reduction efficiency and also the elapsed time gain using GPU acceleration on multi-rate SPIIR filtering. The setup of the performance test is explained in detail here. We simulate two sets of data using Advanced LIGO noise for the twin LIGO detectors, respectively. The duration of each data set is 1000 seconds. These data sets are injected coherently into a simulated binary neutron star coalescence GW signal. To prepare the templates used by our pipeline for the search, we generate two template banks for each detector, which covers the parameters of the injection. Each bank consists of 1024 geometric templates. We generate SPIIR filters for each bank and divided the filters into four groups corresponding to four designated sub-rates. We insert filters with zero coefficients into each group so that the number of filters of each group in a bank will be a power of 2. While it is not necessary to do the filter insertion for the search, the insertion here is for performance testing purpose. The number of SPIIR filters of each group in each bank is shown in table 5.

Table 5: Number of filters in each frequency band of a bank and the data sampling rates required by Nyquist theorem.
Table 6: Performance of the low-latency SPIIR pipeline merely using CPU power (CPU pipeline) and with GPU acceleration of multi-rate SPIIR filtering (GPU pipeline).

Tab. 6 shows the CPU resource reduction and elapsed time gain by the GPU-accelerated SPIIR pipeline. It achieves -fold improvement on CPU resource reduction. This is close to our expectation shown earlier. Besides, we have significantly reduce the running time of the pipeline by 11-fold. The difference of the SNRs between the GPU pipeline and the CPU pipeline on the injection event is within for the 10 test runs, and the injection is successfully recovered by the GPU pipeline.

6 Conclusions and Future Work

Low-latency and real-time detections of gravitational-wave (GW) signals are gaining priority for their potential to enable prompt electromagnetic follow-up observations. The low-latency SPIIR detection pipeline is a CBC detection pipeline with that latencies of tens of seconds.

In this paper, we develop the GPU acceleration for the main computational component of the SPIIR pipeline — multi-rate SPIIR filtering. We first improve the GPU optimization of the filtering part, by employing a new kind of GPU thread configuration and exploiting new memory features of Maxwell GPUs. This provides a notable to 10 fold speedup over our previous acceleration on Fermi GPUs. We implement GPU acceleration of resampling and upstream summation parts for the multi-rate filtering procedure. Our tests show the multi-rate filtering has been accelerated over 100 fold given different filtering scenarios. This leads to a 21-fold CPU resource reduction for the entire pipeline and a 11-fold reduction on elapsed time.

We would like to thank Maurice H.P.M. van Putten for discussion and comments on the details of the paper. This research is supported in part by National Natural Science Foundation of China (Grant No. 61440057, 61272087, 61363019 and 61073008), Beijing Natural Science Foundation (Grant No. 4082016 and 4122039), the Sci-Tech Interdisciplinary Innovation and Cooperation Team Program of the Chinese Academy of Sciences, the Specialized Research Fund for State Key Laboratories, and the Australian Research Council Discovery Grants and Future Fellowship programs. QC gratefully acknowledges the support of an International Postgraduate Research Scholarship funded by the Australian government. We thank Andre Fletcher for proofreading a draft of this manuscript.



  1. Contributed equally to this work with Xiangyu Guo, Corresponding Author’s Email:
  2. Corresponding Author’s Email:
  3. Corresponding Author’s Email:
  4. gstreamer library:
  5. gstlal library:
  6. lal library:


  1. B. P. Abbott, R. Abbott, T. D. Abbott, et al. Observation of gravitational waves from a binary black hole merger. Phys. Rev. Lett., 116:061102, Feb 2016.
  2. Bangalore Suryanarayana Sathyaprakash and Bernard F Schutz. Physics, astrophysics and cosmology with gravitational waves. Living Reviews in Relativity, 12(2):18–19, 2009.
  3. Maurice HPM van Putten, Gyeong Min Lee, Massimo Della Valle, Lorenzo Amati, and Amir Levinson. On the origin of short grbs with extended emission and long grbs without associated sn. Monthly Notices of the Royal Astronomical Society: Letters, 444(1):L58–L62, 2014.
  4. Q. Chu, E. J. Howell, A. Rowlinson, et al. Capturing the electromagnetic counterparts of binary neutron star mergers through low-latency gravitational wave triggers. MNRAS , 459:121–139, June 2016.
  5. B. J. Owen and B. S. Sathyaprakash. Matched filtering of gravitational waves from inspiraling compact binaries: Computational cost and template placement. PRD , 60(2):022002–+, July 1999.
  6. B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton. FINDCHIRP: an algorithm for detection of gravitational waves from inspiraling compact binaries. ArXiv General Relativity and Quantum Cosmology e-prints arXiv:gr-qc/0509116, September 2005.
  7. D. Buskulic, Virgo Collaboration, and LIGO Scientific Collaboration. Very low latency search pipeline for low mass compact binary coalescences in the LIGO S6 and Virgo VSR2 data. Classical and Quantum Gravity, 27(19):194013, October 2010.
  8. K. Cannon, R. Cariou, A. Chapman, et al. Toward Early-Warning Detection of Gravitational Waves from Compact Binary Coalescence. ArXiv e-prints arXiv:1107.2665, July 2011.
  9. J. Luan, S. Hooper, L. Wen, and Y. Chen. Towards low-latency real-time detection of gravitational waves from compact binary coalescences in the era of advanced detectors. ArXiv e-prints 1108.3174, August 2011.
  10. S. Hooper, L. Wen, D. Blair, et al. Low-Latency Detection of Gravitational Waves. In American Institute of Physics Conference Series, volume 1246 of American Institute of Physics Conference Series, pages 211–214, June 2010.
  11. S. Hooper, S. K. Chung, J. Luan, et al. Summed parallel infinite impulse response filters for low-latency detection of chirping gravitational waves. PRD , 86(2):024012, July 2012.
  12. Yuan Liu, Zhihui Du, Shin Kee Chung, et al. Gpu-accelerated low-latency real-time searches for gravitational waves from compact binary coalescence. Classical and Quantum Gravity, 29(23):235018, 2012.
  13. H. S. Black. Modulation Theory. New York, Van Nostrand, 1953.
  14. CUDA Nvidia. Nvidia cuda c programming guide. NVIDIA Corporation, 120, 2011.
  15. Khronos OpenCL Working Group et al. The opencl specification. version, 1(29):8, 2008.
  16. Kate Gregory and Ade Miller. C++ AMP: Accelerated Massive Parallelism with Microsoft® Visual C++®. ” O’Reilly Media, Inc.”, 2012.
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.