An Investigation on Numerical Bugs in GPU Programs Towards Automated Bug Detection

CCS Concepts: • Software and its engineering → Software testing and debugging; Empirical software validation.

Additional Key Words and Phrases: CUDA programming, fault localization, numerical bugs, floating-point errors, bug detection

ACM Reference Format:

Ravishka Rathnasuriya, Nidhi Majoju, Zihe Song, and Wei Yang. 2025. An Investigation on Numerical Bugs in GPU Programs Towards Automated Bug Detection. Proc. ACM Softw. Eng. 2, ISSTA, Article ISSTA073 (July 2025), 24 pages. https://doi.org/10.1145/3728950

Introduction

Numerical applications heavily depend on floating-point arithmetic, which approximates real numbers across varying precision levels, including half, single, and double precision [7]. GPUs are engineered to efficiently execute floating-point operations, offering substantial improvements in computational performance. Notably, NVIDIA GPUs [52] have progressively advanced with each hardware generation. Early NVIDIA architectures supported only single precision, but modern architectures now fully implement IEEE standards for half, single, and double precision, along with

Authors’ Contact Information: Ravishka Rathnasuriya, University of Texas at Dallas, Richardson, USA, rsr200002@ utdallas.edu; Nidhi Majoju, University of Texas at Dallas, Richardson, USA, nidhi.majoju@ utdallas.edu; Zihe Song, University of Texas at Dallas, Richardson, USA, zihe.song@ utdallas.edu; Wei Yang, University of Texas at Dallas, Richardson, USA, wei.yang@ utdallas.edu. Fused-Multiply-Add (FMA) operations introduced in the Fermi [36] and Kepler [37] architectures. These advancements have introduced distinct characteristics in GPU computations compared to CPU computations, including unique rounding mechanisms, specialized mathematical libraries, and differing numerical capabilities (e.g., compute capability versions 2.0, 3.0, 6.x).

There are numerous studies on non-numerical bugs in GPU programs and numerical bugs in CPU programs, but none specifically address GPU numerical bugs (GPU-NBs). Existing research in the GPU domain has primarily focused on detecting synchronization bugs [4, 6, 9, 30, 76], memory bugs [6, 23, 30, 58], and compiler bugs [24]. However, GPU-NBs represent a distinct category arising from the finite approximation of real numbers in computing systems. These bugs can occur due to violations of mathematical properties, numerical operator implementation errors, insufficient data type support, or floating-point representation errors such as round-off and truncation errors. While both numerical and non-numerical bugs can produce similar symptoms, their causes differ. For instance, synchronization bugs result from misuse of barrier APIs [76] causing thread misbehavior, whereas numerical bugs stem from floating-point representation errors.

Despite the substantial body of research on CPU numerical bugs, a dedicated study of GPU-NBs is necessary due to critical distinctions in their manifestation, detection, debugging and repair [3, 8, 11, 60, 71, 72, 78, 80, 82, 83, 84]. First, GPU and CPU architectures differ fundamentally [74], affecting how numerical bugs emerge and are triggered. For instance, unlike CPUs, GPUs lack automated fault tolerance mechanisms such as exception flags and exception status registers [28, 50], which forces developers to use manual techniques (e.g., insert conditional statements) to identify and diagnose errors. Second, CPUs often perform floating-point calculations using an 80-bit “extended” precision mode, whereas GPUs typically use 32-bit and 64-bit for single and double precision operations, respectively [10]. This disparity in precision impacts how rounding-off errors develop, making testing approaches for CPU numerical bugs ineffective for GPUs, which round more frequently [75] and often lack the hardware-level rounding controls [73] available in CPUs. Third, CPUs typically process operations sequentially, facilitating precise tracking of error accumulation, whereas GPUs run thousands of threads concurrently, leading to complex error patterns that require different debugging and detection approaches. Due to these differences, CPU numerical bug detection strategies [3, 8, 11, 34, 57, 60, 71, 72, 80, 82, 83, 84] cannot be directly applied or work effectively in GPU-NB detection. Additionally, CPU strategies for numerical bug fixes often involve simple precision adjustments, but GPU fixes demand custom kernel strategies and precision tuning to maintain performance under parallel execution. Lastly, the proprietary nature of accelerated libraries used in GPU programs [38, 39, 40, 41, 42, 43, 44, 45] restricts access to source code, preventing the use of tools like Clang/LLVM for instrumentation [28]. These factors underscore the need for a focused study on GPU-NBs, examining their symptoms, root causes, and fixing strategies. Such research can improve developers’ understanding, and guide future efforts in automated bug detection for GPUs.

Numerical bugs in GPU programs and deep learning (DL) applications originate from different sources. DL numerical issues arise from flaws in the design and logic of algorithms, such as unstable optimization methods or improper weight initialization strategies, affecting prediction accuracy at the application level [25, 53, 56, 71, 77, 82]. In contrast, GPU-NBs stem from how data types and numerical operations are handled during execution, leading to errors at the hardware level. Rather than analyzing entire applications, our focus is on the numerical correctness of core GPU computations and library implementations.

DL-based detection approaches rely on whole-application testing methods, such as metamorphic and differential testing [71, 77], which perturb models and compare outputs to detect inconsistencies. However, these methods fail to isolate hardware-level numerical errors, as they would require exhaustive exploration of vast model and parameter spaces for whole-application testing. The absence of hardware-level detection techniques leaves a critical gap, as GPU-NBs introduce hidden computation errors without immediate failure signals, making them undetectable through DL-based detection methods.

The purpose of our study is to investigate GPU-NBs, focusing on common issues and solutions related to floating-point operations. We analyzed bug reports from both NVIDIA’s open-source and closed-source numerical libraries, manually examining 906 GPU bug samples, of which 397 were identified as numerical issues. Based on these findings, we categorized numerical bugs into four main types: precision bugs, special-value bugs, device-specific bugs, and data type and operator bugs. We discuss the characteristics of each category, including their symptoms, test input patterns, oracles that trigger the bugs, and the strategies used to fix these numerical issues in GPU programs.

Our study offers findings and insights that pave the way for future researchers to develop a systematic framework for addressing GPU-NBs. The framework begins with bug detection (§ 5), encompassing test input generation and the selection of test oracles. Based on symptoms identified during detection, the framework advances to debugging (§ 4) by mapping symptoms to potential root causes (Figure 2). Once the root causes are identified, the framework applies corresponding fixing strategies (Figure 2) to repair the program (§ 6). We propose a preliminary tool, GPU-NBDetect, for detecting GPU-NBs (§ 5.3). GPU-NBDetect employs seven mutation rules to generate floating-point inputs and detects numerical bugs across six categories. GPU-NBDetect identified 226 bugs across 186 mathematical functions in four libraries, with 60 confirmed by developers, outperforming four existing bug detectors and two input generation methods.

Our main contributions of this study are:

• Summarizing patterns of the symptoms, the root causes, and the fixing strategies of numerical bugs in GPU programs.

• Providing guidance for creating reliable and stable GPU program test inputs and test oracles for detecting numerical bugs.

• Outlining common fix strategies for numerical issues, helping developers reproduce and resolve bugs efficiently, thus conserving development resources.

• Offering insights into the automated detection, debugging, and fixing of GPU-NBs.

• Introducing GPU-NBDetect, a preliminary tool featuring input generation and test oracle selection to detect GPU-NBs.

Background and Motivation

Floating-point numbers approximate real values based on the IEEE-754 standard [7]. GPUs support multiple formats: half, single, double, tensor and brain floating point [17], making format selection critical for ensuring computational accuracy in applications ranging from scientific computing to DL.

While GPUs excel at parallel computation, they are prone to numerical inaccuracies, particularly in floating-point arithmetic. Figure 1 highlights discrepancies in reduction sum operation results between CPU and GPU computations for values like 0.0008 and 0.4999. For example, the CPU computes 799.969, while CUDA yields 794.618; similarly, 499999.500 on the CPU becomes 499884.438 on the GPU. These deviations stem from cumulative rounding errors [17, 74, 75] and non-associative nature of floating-point operations, reflecting the precision limitations of GPU computations.

These examples highlight a broader issue: as computation complexity increases, GPU-NBs can cause substantial deviations from expected results, which is critical for accuracy-dependent applications. Despite the manual developer efforts to adjust precision or error-checking mechanisms, the root causes of GPU-specific inaccuracies remain insufficiently explored. This study offers a comprehensive analysis of GPU-NBs, examining real-world cases to identify root causes, symptoms, input patterns, bug-triggering oracles, and effective mitigation strategies.

Methodology

In this section, we describe our data collection process (§ 3.1) and discuss the methodology for filtering and labeling GPU-NBs (§ 3.2 and § 3.3).

Data Collection

Building a comprehensive dataset of GPU-NBs requires systematic identification of relevant issues across libraries and projects. Since GPU libraries lack standardized popularity metrics like power scoring [19], we used GitHub as the primary source, collecting bug reports from January 2020 to June 2024 through a combination of automated queries and manual review to ensure relevance.

Our methodology distinguishes between closed-source and open-source GPU libraries. For closed-source libraries, such as cuFFT [41], cuBLAS [40], Math [39], cuRAND [42], cuSOLVER [43], cuSPARSE [44], and cuTENSOR [45], we collected bug samples from dependent projects due to limited direct access. In contrast, open-source libraries like CuPy [61], Numba [68], MatX [46], and Ginkgo [12] allowed direct analysis of the libraries and their associated projects. This approach enabled the construction of a diverse GPU-NB dataset spanning domains from scientific computing to DL.

Issue Selection

Identifying GPU-NBs across GitHub is time and labor-intensive, therefore, we adopted a systematic keyword-based filtering approach tailored to both closed- and open-source libraries. This section outlines the criteria and methods used to extract relevant bugs and eliminate invalid reports.

Selection GPU-NBs related to closed-source libraries. Using the GitHub search API [13], we initially collected 118,438 reports with broad keywords such as “GPU” and “GPU programming.” To focus on numerical bugs, we applied a refined set of keywords based on methodologies from previous studies on CPU numerical bugs [8]. These keywords included terms such as “nan”, “overflow”, “underflow”, “infinity”, “precision”, “unstable”, “roundoff”, “truncation”, “cancellation”, “accuracy”, “numeric”, “floating”, “fp”, “cufft”, “cublas”, “curand”, “cuda math library”, “cusolver”, and “cusparse”. This filtering process narrowed our sample to 10,380 issues.

Removing Invalid Issues: To ensure relevance and reduce noise, we applied further criteria: (1) confirmed bugs (closed status), (2) consists of a documented fix, and (3) classified as a numerical bug (excluding memory and synchronization issues). These steps narrowed the sample size to 624.

Selection GPU-NBs related to open-source libraries. For open-source libraries, we expanded the search to include both the libraries and their dependent projects. Using terms like “CUDA” and “numerical”, we initially identified 76 repositories. We then combined “CUDA” with key numerical libraries from Franco et al. [8], yielding additional repositories for “NumPy”, “SciPy”, “LAPACK”, “GSL”, and “Elemental”. After removing duplicates and irrelevant projects, we finalized a set of open-source libraries and projects, including CuPy [61], Numba [68], MatX [46], and Ginkgo [12].

Removing invalid issues: We refined the dataset by retaining only closed bug reports, as open reports may reflect unresolved or unverified bugs. This final filtering yielded 282 validated GPU-NB cases.

GPU-NB Labeling

After merging both datasets, removing duplicates, and excluding misclassified bug reports, we obtained 397 unique GPU-NB cases (229 from closed-source and 168 from open-source libraries). These were categorized by root cause using a three-step process.

First, two authors independently reviewed each bug report by analyzing descriptions, comments, and associated commits. While we initially applied Franco et al.’s CPU-focused framework [8], it lacked specificity for GPU-related issues. Notably, the CPU categories for accuracy and convergence overlapped, as both addressed round-off and truncation errors. To resolve this, we merged them into a broader precision bugs category, which we expanded to include arithmetic, approximation, and casting bugs unique to GPUs. Additionally, we introduced two new categories: data type and operator bugs and device-specific bugs to capture GPU-architecture-specific issues. This refinement resulted in seven distinct categories, with only special value bugs overlapping with CPU studies.

Second, we ensured consistency by comparing initial categorizations and merging similar categories as needed. A joint review was conducted to resolve 25 inconsistently labeled reports.

Finally, we mapped specific symptoms to each root cause, enabling a fine-grained categorization of GPU-NB manifestations. Unlike the broader symptom types used in CPU studies, our detailed classification captures GPU-specific behaviors and reveals strong one-to-one correspondences between symptoms, root causes, and fixes. We further analyzed common repair strategies, noting a reliance on manual fixes in GPU development, in contrast to the automation observed in CPU-focused work. The complete mapping is illustrated in Figure 2, and a summary of the analyzed GitHub reports is available on our project website [18].

Causes of Numerical Bugs

We identified 397 numerical bugs, which we categorized into four primary groups: precision, data type and operator, special value, and device-specific bugs. Table 2 presents the distribution of each bug type across these categories.

Root Cause	Sub categories	Closed-Source	Open-Source	Total
Precision Bugs	Arithmetic Operation Bugs	20	11	31
Precision Bugs	Rounding/Approximation Bugs	25	9	34
Precision Bugs	Casting Bugs	23	9	32
Data Type and Operator Bugs	Data Type Bugs	84	59	143
Data Type and Operator Bugs	Operator Bugs	37	57	94
Special Value Bugs	Special Value Bugs	33	15	48
Device-specific Bugs	Device-specific Bugs	7	8	15

Precision Bugs

Precision bugs in GPU programs arise from limitations in representing floating-point numbers accurately, leading to computations with insufficient or imprecise numerical values. Approximately 24.4% of the identified numerical bugs fall into this precision category (I-1 in Table 1). These bugs fall into one of three subcategories:

	Findings about numerical bug causes	Implications on unique characteristics of GPU-NBs
1	Among the 397 observed GPU program tests, 97 bugs were identified as precision bugs, including overflow/underflow, rounding-off, and casting errors. (Table 2)	Precision bugs in GPUs occur because, unlike CPUs with extended precision (80-bit mode) and automatic error detection, GPUs lack mechanisms to control precision, causing errors to propagate unnoticed.
2	One-third of precision issues are due to rounding errors (Table 2).	Rounding errors arise because GPUs lack programmable rounding modes and extended precision support. Unlike CPUs, which allow flexible precision adjustments, GPUs use fixed rounding mechanisms, leading to cumulative precision loss, particularly in parallel computations.
3	One-third of precision issues are due to casting errors (Table 2). Studies on CPU programs have not identified casting errors.	Casting errors occur because, unlike CPUs, GPUs substitute default values for unsupported data types rather than raising exceptions, allowing errors to go undetected.
4	Among the 397 observed GPU program tests, 239 (59.6%) are caused by data type and operator issues unique to GPU, particularly with CUDA numerical operations (Table 2).	GPU software bugs stem from unsupported data types and unimplemented operations in CUDA, as limited transparency and testing in GPU libraries leave compatibility issues undetected, unlike fully verified CPU libraries.
5	Special value bugs account for 12.1% of the observed GPU program tests (Table 2).	GPU bugs arise from a lack of automatic NaN and INF detection, causing silent accumulation of rounding errors. Unlike CPUs with built-in exception alerts, GPUs rely on manual checks to handle these inconsistencies.
6	Environmental issues account for 3.7% of numerical problems in GPU programs and are unique to GPU (Table 2).	GPU bugs stem from inconsistent precision support across device generations. Designed for speed over accuracy, GPUs lack uniform numerical support, forcing developers to manually manage precision inconsistencies.
	Findings about numerical bug manifestation	Implications on bug detection
7	Developers use floating-point numbers, integers, out-of-range values, and NaNs as inputs to GPU programs (§ 5.1).	Automated, feedback-driven fuzzing improves GPU bug detection by dynamically generating test inputs that expose numerical bugs, addressing GPU’s sensitivity to data types and rounding behavior.
8	All GPU programs can be tested against their corresponding CPU versions (Figure 10).	Automated cross-platform testing enhances GPU bug detection by using CPU benchmarks to reveal GPU-specific discrepancies, leveraging the CPU’s robust precision control.
9	Of the 397 tests, 179 (45.1%) have the capability to log program results (Table 3 and Table 4).	Automated logging aids GPU bug detection by capturing results and comparing them with CPU outputs, effectively identifying silent rounding and casting issues.
10	Of the 397 tests, 37 (9.3%) can profile program performance in terms of execution (Table 3 and Table 4).	Automated performance profiling can enhance GPU bug detection by flagging bottlenecks linked to precision-speed trade-offs, which are less pronounced in CPU.
	Findings about numerical bug symptoms	Implications on debugging
11	The observed GPU program tests exhibit three primary symptoms: error messages (45.6%), wrong results (45.1%), and performance bottlenecks (9.3%) (Figure 2 and Table 3).	Debugging GPU-NBs reveals symptoms closely tied to root causes (Figure 2). Unlike CPUs, GPUs experience more severe performance bottlenecks from numerical bugs, often leading to significant performance losses (Table 4).
12	Error messages are prominent in GPU programs, occurring in 181 out of 397 cases (Table 3).	Unique error messages in GPUs such as overflow/underflow, invalid casting, unsupported data types, operators, and undefined computations can directly reveal root causes (Fig. 2). Unlike CPUs with fault-tolerant mechanisms like exception detection, GPUs require developers to manually insert such checks for effective debugging.
13	Producing wrong results helps identify numerical bugs in GPU programs. Our study found that 179 out of 397 issues (45.1%) were confirmed by comparing program results (Table 3).	Debugging wrong results in GPUs involves comparing outputs across scenarios, such as varying precision levels, casting operations, data types, and checking for NaN values (Figure 2).
14	Performance bottlenecks help detect numerical bugs even when programs produce correct results, occurring in 9.3% of cases (Table 3).	Debugging performance bottlenecks in GPUs involves comparing cast operations, data types, and CPU execution times, with visual profilers aiding analysis of time and memory usage (Figure 2).
	Findings about numerical bug fixing	Implications on program repair
15	Various precision issues are addressed through different methods: 27.8% with checkers, 39.3% by changing precision, and 32.9% using upcasting/downcasting strategies (Table 7 and Figure 2).	Automated repair approaches for precision issues could save time by adjusting numerical precision without relying on manual, trial-and-error methods, which are typically time-consuming for users.
16	Four out of 57 (7.1%) operator issues in third-party libraries are resolved using custom kernel strategies, which are unique to GPU programs (Table 7 and Figure 2).	CuPy’s custom kernel strategy assists users in defining numerical operators within GPU programs.
17	Data type issues in GPU programs can only be resolved by adding data type support or using a compatible data type (Table 7 and Figure 2).	Automated repair could address data type issues by adding support for required types or selecting compatible data types, reducing the need for manual adjustments.

(1) Arithmetic Operation Bugs: Arithmetic operation bugs arise when computations exceed or fall below the supported numerical range, resulting in overflow or underflow. Unlike CPUs, which use exception status registers (e.g., fenv.h in C) to detect such conditions, GPUs lack built-in mechanisms for numerical exception handling. As a result, computations may silently yield values like +inf, -inf, or 0.0 [28] unless developers manually implement detection logic.

These bugs are particularly prominent in GPU programs due to the frequent use of mixed-precision computation [14, 16] and integer quantization [15]. For example, values stored in FP32 may be processed in FP16, increasing the risk of range violations. Additionally, implicit conversions from integers to floating-point types can introduce further precision loss and memory overhead. As shown in Figure 2, arithmetic operation bugs in GPU code often lead to specific error messages, such as developer-inserted overflow or underflow flags.

An example from the CuPy project [65] illustrates this problem in Figure 3. During the convolve2d operation, an integer overflow occurs because the operation is executed with incompatible data types. When applied to two integer arrays, $a$ and $b$, with different precisions (uint64 and uint8, as shown in lines 4 and 5 in Figure 3a), the SciPy library correctly handles the operation in the higher precision data type (uint64), avoiding overflow (line 6). In contrast, CuPy defaults to the lower precision type, uint8 (line 15 in Figure 3b), resulting in integer overflow since uint64 values cannot be accurately represented in uint8.

(2) Rounding/Approximation Bugs: Bugs in this category occur when real numbers are approximated within finite precision limits, causing deviations from expected values. While both CPU and GPU programs are affected, the impact is more pronounced on GPUs due to their limited precision and lack of extended precision support. Precision loss is further exacerbated by operation ordering effects such as subtractive cancellation [35].

CPUs often mitigate such errors using 80-bit extended precision, whereas GPUs are restricted to 32-bit or 64-bit precision and lack configurable rounding modes [17, 74] (I-2 in Table 1). Instead, developers must rely on compiler intrinsics to control rounding behavior. As shown in Figure 2, rounding bugs in GPUs can lead to significant discrepancies between actual and expected outputs.

An example of this can be seen in the CuPy project [64] in Figure 4, where the map_coordinate operation with order=none and order=0 produces incorrect outputs compared to SciPy. In this case, the input array point in line 2 contains a value of 0.5, leading to a rounding error in CuPy’s precision handling. CuPy uses the rint() function, which rounds to the nearest integer, whereas SciPy employs floor(double), which rounds down to the nearest double value. This difference in rounding methods leads to output inconsistencies between the two libraries.

(3) Casting Bugs: Casting bugs occur when values are converted between data types, resulting in numerical inconsistencies. While relatively rare in CPU numerical libraries [8], casting issues are more prevalent in GPU programs, where such operations often produce inconsistent or erroneous outcomes (I-3 in Table 1).

Casting between high and low precision data types introduces risks of precision loss, which should ideally remain within acceptable bounds (e.g., within 0.1). Unlike CPUs that leverage register flags to detect casting errors, GPUs lack exception-handling mechanisms. As a result, when a value is cast to an unsupported type, GPUs often substitute default values, allowing errors to propagate silently and undetected.

As shown in Figure 2, casting issues in GPUs can be detected by (1) examining discrepancies between numerical values before and after casting, and (2) identifying performance bottlenecks when converting between precision types.

An example of a casting issue is shown in the cuDF project [49], illustrated in Figure 5. In this program, the unary.cast operation between Decimal64 and int32 yields erroneous results (lines 6 and 9). The unary.cast operator lacks an exception detection mechanism, failing to verify whether the source value can be accurately represented in the target data type, int32_t. Consequently, casting -72462.12 from Decimal64 to int32 produces an inaccurate default value of 13437, while casting -48938 from int32 to Decimal64 results in -5988.32704.

Furthermore, a cuDF example [47], shown in Figure 6a and 6b illustrates how data type choice affects casting performance. Casting from Decimal64 to Decimal128 required 4134 ms, while casting between identical Decimal128 types took only 1598 ms, indicating a 2.5-fold difference in execution time. In these figures, grey sections indicate the casting operation duration, and green sections represent CUDA API activities.

Data Type and Operator Bugs

Data types [5] specify how values are represented, while operators [70] define how computations are performed on them. Bugs related to data types and operators account for 59.6% of numerical issues in GPU programs (I-4 in Table 1). These bugs typically arise from insufficient support for certain data types or operations, resulting in incorrect behavior or performance degradation.

Data type bugs arise when unsupported or incorrect types are used in computations, often leading to invalid casting or unsupported data type errors (Figure 2). For example, FP16 operations may perform worse than FP32 under certain conditions, causing unexpected performance bottlenecks. Operator bugs result from missing or incorrectly implemented functions, leading to unsupported operations or erroneous outputs. These two categories are often interrelated, as data type issues can be misinterpreted as operator failures.

While GPU libraries often adopt CPU-style implementations for numerical operators, they differ significantly in documentation quality and transparency. Many are closed-source, preventing developers from verifying supported data types or internal behaviors. Consequently, operators that work correctly on CPUs may fail on GPUs due to inadequate type support. This opacity complicates debugging, as users cannot inspect or modify internal implementations. Unlike CPU libraries, which primarily address typographical issues [8], GPU-NBs more frequently stem from unsupported data types and operator incompatibilities.

In the cuDF project [48], an example of a data type bug occurs when creating a series of float16 CuPy objects, causing a crash with an “invalid data type” error in Figure 7. This issue arises because the __cuda_array_interface__ expects np.dtype objects, but the program incorrectly uses cudf.dtype, leading to a failure.

An operator-related issue is observed in the CuPy project [62] in Figure 9, where the program fails to pass the cupy.mean function into the generic_filter() operation, resulting in a “bad type error” message. However, the equivalent np.mean function in NumPy passes without issue (line 3), highlighting a compatibility problem with CuPy’s handling of certain operators.

Special Value Bugs

Special value bugs arise when computations produce or encounter NaN or (+/-)INF values, accounting for 12.1% of the dataset (I-5 in Table 1). These bugs often occur when numerical operations yield undefined results, causing error messages or incorrect outputs due to the presence of NaN or INF values (Figure 1).

A primary cause of special value bugs in GPUs is the lack of built-in exception handling for NaN and INF values. Unlike CPUs, which include exception detection for these cases, GPUs require developers to manually implement checks for special values. Without automated detection, errors involving NaN or INF can propagate undetected, resulting in incorrect computations.

Differences in NaN handling between CUDA and CPU-based libraries have also been observed [55]. For example, NVIDIA GPUs use little-endian format exclusively, whereas CPU-based libraries like NumPy support both little- and big-endian formats. CuPy, however, limits support to little-endian, which can lead to incorrect results when processing big-endian data [66].

An example from the PyTorch project [55] illustrates a special value bug (Figure 8). NaN values appear in a.grad computations with FP16 values (line 5) when a small value is passed to the torch.tensor() function. This results in INF values that propagate and produce NaN outputs.

Device-specific Bugs

Device-specific bugs in GPU programs arise from variations in numerical support across different GPU models. As shown in Table 2, 15 out of the 397 identified bugs are attributed to device-specific factors (I-6 in Table 1). These bugs often result in performance delays compared to CPU execution time when using GPUs with limited numerical support (Figure 2).

While CPU architectures generally maintain backward compatibility and consistent numerical support across generations, GPU support can vary, leading to device-specific issues [16]. For example, older GPUs lack support for mixed-precision computing, while newer models, such as the RTX series, V100, and A100, achieve substantial performance improvements through mixed precision.

In the CuPy project [63], performance discrepancies in double precision computations illustrate this issue. Tests reveal that matrix multiplication with NumPy on a CPU is faster than the CuPy implementation on a GPU. Specifically, the CPU execution time with NumPy was 0.7268 seconds, while the GPU execution time with CuPy was 1.0382 seconds. The test was conducted on a GeForce GTX 760 GPU, which is known for its low double precision performance.

Mapping Symptoms to Root Causes

This section summarizes methods for identifying the root causes of GPU-NBs based on their symptoms. Our study identifies three symptom types resulting from GPU-NBs: incorrect results, error messages, and performance bottlenecks, as depicted in Figure 2. Moreover, we extend the classification with detailed subcategories for each symptom type. Table 3 presents the distribution of these symptoms across identified GPU-NBs.

Root Cause	Sub categories	Closed-source	Open-source	Total
Error Messages	Overflow/Underflow	21	15	36
	Casting	3	8	11
	Data Type	57	33	90
	Operator	19	18	37
	Undefined Computation	4	3	7
Wrong Results	Before and After Casting	17	3	20
	Different Data Types	9	17	26
	Different Precision Range	16	34	50
	NaN Outputs	30	13	43
	Different Amount of Precision	27	13	40
Performance Bottlenecks	Cast Operations	2	0	2
	Data Types	17	4	21
	CPU Programs	7	7	14

Unique error messages directly indicate the root cause of a GPU-NB, particularly when developers place manual checkers in GPU code due to the lack of exception detection capabilities. For example, overflow or underflow errors suggest that a computed value has exceeded the supported precision range, indicating an arithmetic operation bug. Errors such as “invalid casting” occur when a numerical value is cast to an unsupported data type, while “unsupported data type” or “unsupported operator” errors signal attempts to use data types or operations not supported by the API. Special value bugs are typically flagged by errors involving undefined values like NaN or INF (F-12 and I-12 in Table 1).

Comparing GPU and CPU outputs for the same program provides further insight into GPU-NBs. Manually adapting GPU programs to corresponding CPU implementations (e.g., replacing CuPy with NumPy or SciPy) helps detect discrepancies. Casting bugs are confirmed when results differ before and after casting, while precision discrepancies point to rounding bugs. Varying results across different data types for the same values suggest a data type bugs, and differences across precision ranges for specific operators indicate an operator issue. NaN outputs consistently identify special value bugs (F-13 and I-13 in Table 1).

For programs yielding correct results, performance measurement can also reveal numerical bugs. Performance bottlenecks often indicate underlying issues; for instance, casting operations can introduce inefficiencies, with lower precision types like FP16 sometimes performing worse than FP32. Programs that show poorer performance on GPUs than CPUs may signal device-specific bugs (F-14 and I-14 in Table 1), as detailed in § 4.

Bug Detection Strategies

Our study highlights numerical bug detection strategies, as summarized in Figure 10. Based on findings from Table 1, we identify two primary components in detecting GPU-NBs: test input generation and test oracles. The goal of test input generation is to create inputs that can trigger numerical bugs. We explored two methods: Predefined General Inputs and Feedback-Driven Customized Input. Once inputs are generated, the second stage involves applying test oracles to identify GPU-NBs.

Test oracles use output comparisons to detect bugs, with some strategies requiring single test runs (e.g., Strategies 1-5 and 9 in Table 4), while others rely on multiple runs.

	Program Iterations (input generations, programs)	Manifestation Strategy	Assertion	Root Cause
1	Single input, single GPU iteration	Logging result A of GPU version	if A > MAX(FP) or A < -MAX(FP)	Overflow bug
2	Single input, single GPU iteration	Logging result A of GPU version	if A < MIN(FP) or A > -MIN(FP)	Underflow bug
3	Single input, single GPU iteration	Logging result A of GPU version	if A == NaN	Special value bug
4	Single input, single CPU and GPU iteration	Logging precision in result A of CPU version and result B of GPU version	abs(A-B) > 1E-7	Rounding bug
5	Single input, single GPU iteration	Logging value A before casting and value B after casting	abs(B-A) > 0.9	Casting bug
6	Multiple inputs, multiple GPU iterations	Logging result A of first run and result B of another run with different data types for the same value	abs(A-B) > 0.01	Data type bug
7	Multiple inputs, multiple GPU iterations	Profiling execution time A of cast operation and time B of cast operation with higher data type	lower data type A slower than B (avg 2.5X slower)	Casting error bug
8	Multiple inputs, multiple GPU iterations	Profiling overall time A of one data type and time B of higher data type	A is slower than B (nearly 3x slower)	Data type bug
9	Single input, single CPU and GPU iterations	Profiling time A of CPU version and time B of GPU version	A is faster than B	Device-specific bug

Test Input Generation

Our findings indicate that many developers rely on random input modifications to trigger numerical bugs in GPU programs. Developers start with predefined inputs and apply random perturbations, such as altering precision, to generate new inputs. This ad hoc approach can be systematically automated using fuzzing techniques, which offer a structured means of input variation to improve bug detection (F-7 and I-7 in Table 1).

5.1.1 Predefined General Inputs. These inputs serve as the initial test values used to run the program and are categorized into four types: out-of-range numbers, integers, floating-point numbers, and NaN. Out-of-range numbers fall outside the representable range of a given data type, such as the range for 32-bit floating-point numbers, which spans from 1.17549435E-38 to 3.40282347E+38 [22]. NaN values represent undefined or unrepresentable results, often arising from invalid computations or missing inputs [1]. Floating-point numbers of various precision levels, such as half, single, double, and decimal, are used to test for precision-related issues. Additionally, integers serve as test inputs due to CUDA’s implicit conversion of integers to floats, which may introduce precision errors. This categorization aids in understanding how distinct input types impact numerical accuracy in GPU computations.

5.1.2 Feedback Driven Customized Inputs. These inputs leverage mutation techniques to systematically generate inputs that trigger numerical bugs. These inputs are categorized into value-based and type-based mutations. Value-based mutation employs fuzzing techniques [81] to systematically vary numerical values, increasing the likelihood of identifying GPU-NBs. Type-based mutation [81] introduces diverse data types to verify the program’s handling of type-specific issues. Advanced methods like Greybox fuzzing [81] can be employed to apply these mutations, automating the process of discovering new inputs that might trigger hidden bugs.

Test Oracles

Test oracles are essential for verifying the correctness of results in GPU programs. As shown in Figure 10, our study finds that correctness is initially confirmed by comparing GPU results with outputs from corresponding CPU programs using the same predefined inputs. When discrepancies arise, customized inputs or manifestation strategies can be applied to further investigate numerical bugs (F-8 and I-8 in Table 1). These strategies rely on assertion thresholds to differentiate expected numerical variations from valid bugs, and thus require carefully calibrated error bounds to ensure reliable detection.

To define these thresholds, we conducted an empirical evaluation across a diverse set of GPU programs and input types. For each manifestation strategy (Table 4), we measured the distribution of deviations observed between CPU and GPU outputs, or across data type and performance profiles. Thresholds such as 1E-7 for rounding bugs and 0.9 for casting bugs were selected to capture meaningful numerical divergence while minimizing false positives. Similarly, thresholds like 0.01 for data type bugs and 2.5x–3x slowdowns for performance profiling were informed by consistent performance gaps observed across multiple devices and libraries. These values were validated through repeated testing and manual review, ensuring that the selected bounds were robust across architectures and sensitive enough to detect GPU-NBs.

When result discrepancies are observed, error messages offer crucial insight into the root cause, as each numerical bug typically generates a distinctive message. Strategies 1-3 in Table 4 provide guidance on implementing exceptions for detecting overflow, underflow, and special value bugs. If no error message is present, alternative manifestation strategies based on the mutation technique used to generate test inputs can be employed.

For value-based mutation, logging and profiling mechanisms can locate numerical bugs. Strategies 4 and 5 in Table 4 outline approaches for detecting incorrect results. Rounding bugs, for example, can be identified by comparing logged precision levels between CPU and GPU outputs (e.g., Figure 4), while casting bugs are detected when pre- and post-cast values differ (e.g., Figure 5). Profiling can reveal performance bottlenecks, as casting to lower data types may unexpectedly perform slower than higher types (e.g., Figure 6a and Figure 6b). Device-specific bugs become evident when GPU execution time is slower than CPU execution, as seen in CuPy [63]. If all other strategies yield consistent results, users can test with varying precision ranges to detect potential operator issues.

For type-based mutation, logging and profiling are also effective tools (strategies 6 and 8 in Table 4). Programs are executed with multiple data types on the same input values, requiring at least two runs to enable comparative analysis. Data type bugs are identified when program results diverge significantly beyond minor tolerances (e.g., exceeding a difference of 0.01). Profiling execution times across data types can further reveal data type bugs, especially when operations on a lower precision type take longer than those on a higher precision type.

Preliminary Tool

Based on the findings from § 5.1 and § 5.2, we developed GPU-Numerical Bug Detector (GPU-NBDetect) to detect GPU-NBs. GPU-NBDetect integrates automated input generation with differential testing against CPU programs to validate numerical correctness in GPU computations. Its detection logic is guided by the manifestation strategies outlined in Table 4 and is informed by key insights from Table 1, particularly the importance of input generation (F-7, I-7) and cross-platform output comparison (F-8, I-8) for uncovering GPU-NBs.

GPU-NBDetect uses a template-based approach to automatically generate CUDA and C code for testing. The tool dynamically creates the necessary files in Python, requiring users to simply specify the target mathematical APIs. For CPU programs, the templates include predefined operators, input argument handling, and essential libraries such as math.h for numerical operations and exception handling. The GPU templates cover both host and device code, managing memory allocation, data transfer, and kernel execution. Differential testing is performed by executing each compiled program through Python’s subprocess interface.

GPU-NBDetect begins testing with user-defined inputs, initialized by selecting random values within the supported range of each mathematical operator, based on documentation [39]. To efficiently explore the vast input space and uncover bugs, the tool employs seven mutation operators inspired by fuzzing techniques [33]. These mutators generate diverse numerical inputs, particularly near boundary values (e.g., 0.5, 0.499999, 0.50000001), reflecting the common developer practice of manually probing edge cases to reveal numerical issues. The tool employs the following mutators:

• bitflip_s: Flips the sign component of a floating-point number with a 50% probability, switching between 0 and 1.

• bitflip_e: Flips a single bit at a random position within the exponent.

• bitflip_m: Flips a single bit at a random position within the mantissa.

• add_e: Adds a randomly generated integer (between 1 and 32) to the exponent.

• add_m: Adds an integer $2^p - 1$ to the mantissa, where p is randomly generated between 1 and 33.

• sub_e: Subtracts a randomly generated integer (between 1 and 32) from the exponent.

• sub_m: Subtracts an integer $2^p - 1$ from the mantissa, where p is randomly generated between 1 and 33.

• random_bit: Sets a randomly chosen bit in the exponent or mantissa to a random value.

To detect discrepancies, GPU-NBDetect first applies differential testing by comparing GPU outputs with CPU results. Detected bugs are then classified using assertion-based manifestation strategies (Table 4), which target overflow, underflow, rounding, special values, casting, and data type bugs. To ensure robust and accurate classification, multiple error bounds were empirically tested to establish effective thresholds for each strategy.

Minimizing false positives. To reduce false positives, GPU-NBDetect applies input constraints by (1) restricting values to supported floating-point ranges and (2) excluding known bug-triggering inputs identified in official documentation (e.g., CUDA Math Library [39]). Additionally, multiple error bounds were tested to calibrate assertion thresholds for each manifestation strategy (Table 4). Bug samples that failed to consistently meet these thresholds were excluded, improving precision. While this approach ensures accurate bug identification through meaningful result comparisons, profiling-based strategies (strategies 7–9) may still produce occasional false positives due to tool sensitivity and broader input variability.

Evaluation Setup

We conducted experiments on two types of GPU programming environments: C-based and Python-based GPU libraries. For C-based libraries, we utilized CUDA Math [39] and CUDA Thrust [51]. To benchmark the performance of GPU-NBDetect, we compared it against four established C-based GPU tools: XScope [27], which employs Bayesian optimization for input generation; XScope-RS (Random Sampling) [27], which utilizes random sampling techniques; BinFPE [28], which employs the NVBit dynamic binary instrumentation framework to inspect machine registers for exceptions following each computation; and GPU-FPX [32], another NVBit-based tool that detects and analyzes floating-point exceptions on NVIDIA GPUs via binary instrumentation. Since BinFPE [28] and GPU-FPX [32] rely on predefined inputs and do not incorporate advanced input generation techniques, we employed a random input generation strategy akin to random sampling to ensure a fair comparative analysis.

For Python-based GPU libraries, we evaluated CuPy [61] and Numba [68]. Notably, GPU-NBDetect is the first approach to systematically test Python-based CUDA math APIs for numerical bugs. Due to the lack of direct baselines for comparison in Python-based GPU environments, we adopted different input generation strategies. Specifically, we leveraged Bayesian optimization (BO) and random sampling (RS) [27] techniques to evaluate the effectiveness of GPU-NBDetect.

To maintain a focus on mathematical computations, we excluded APIs that do not perform numerical operations (e.g., Min, Max) across all evaluated libraries. We tested 53 functions in both CUDA Math and Thrust, and 40 functions in both CuPy and Numba. Furthermore, we did not include DL-based numerical bug detection approaches [53, 77], as these methods primarily focus on model mutation or hidden-state mutation rather than input mutation. However, XScope-RS employs a strategy similar to the tensor partitioning technique [71], where input components—mantissa, exponent, and sign—are randomly sampled. The key advantage of GPU-NBDetect lies in its robust input generation, mutation, and manifestation strategies, which enable more comprehensive numerical bug detection across a wide range of GPU libraries.

Experiments were conducted on a server with an Intel Xeon E5-26 CPU and eight NVIDIA 1080 Ti GPUs, using CUDA version 11.1. For each mathematical function, GPU-NBDetect generates 100 input samples, with three independent test cycles ensuring result consistency and minimizing the randomness from mutation operators.

Results and Analysis

Table 5 and Table 6 summarize GPU-NBs identified by GPU-NBDetect, which detected a total of 226 bugs across 186 mathematical functions in four libraries, with 60 confirmed by developers. We measure the average time overhead for bug detection in seconds.

Library	Tool	Overflow Bug	Underflow Bug	Rounding Bug	Casting Bug	Data Type Bug	Special Value Bug	Total Bugs	Overhead
CUDA Math	XScope	8	8	-	-	-	12	28	8.7
	XScope-RS	4	5	-	-	-	4	13	2.7
	BinFPE	5	2	-	-	-	5	12	3.2
	GPU-FPX	6	4	-	-	-	5	15	4.2
	GPU-NBDetect	10	14	24	10	14	12	84	1.1
CUDA Thrust	XScope	4	7	-	-	-	6	17	8.4
	XScope-RS	4	4	-	-	-	5	13	3.4
	BinFPE	5	6	-	-	-	8	19	3.9
	GPU-FPX	5	8	-	-	-	9	22	5.0
	GPU-NBDetect	8	16	20	4	8	14	70	1.3

Library	Tool	Overflow bug	Underflow bug	Rounding bug	Casting bug	Data Type bug	Special Value bug	Total bugs	Overhead
CuPy	RS	6	4	2	0	2	4	16	2.1
	BO	7	4	2	0	2	4	19	2.2
	GPU-NBDetect	9	7	4	3	4	8	35	1.7
Numba	RS	6	3	2	1	3	4	19	2.3
	BO	7	5	3	2	3	6	26	2.3
	GPU-NBDetect	9	8	4	3	4	9	37	1.6

Bug confirmation process. We reported the identified bugs on the NVIDIA developer forum, where NVIDIA developers reviewed and verified the findings. Although CUDA documentation (e.g., Math Library [39]) acknowledges that certain inputs may cause bugs, it does not clarify whether these occur within the supported input range. GPU-NBDetect identified these ambiguous cases, prompting developers to create unit tests for validation. Developers confirmed that GPUs lack built-in exception handling and employ different rounding modes compared to C programs. Notably, GPU-NBDetect pinpointed bug-triggering inputs within the valid input range.

GPU-NBDetect reveals six types of numerical bugs through its analysis. GPU-NBDetect identified rounding errors in 52 functions (27.9%), highlighting significant precision issues. Developers remarked that “although the results are within the supported ULP range, these deviations represent errors from the true values.” Notably, rounding bugs were more stable in CuPy and Numba than in C-based GPU libraries, suggesting that Python-based libraries have more robust rounding modes. Overflow, underflow, and special value bugs were also prevalent but challenging to detect due to the absence of floating-point exception handling in GPUs. As developers clarified, “CUDA device code does not handle exceptions, so it’s essential to verify the behavior against C code.” Casting bugs were detected in 14 operations, indicating that GPU casting mechanisms are particularly error-prone across both C-based and Python-based libraries. Furthermore, data type bugs emerged when lower precision values led to less efficient calculations compared to higher precision numbers.

In comparison with C-based CUDA libraries and other tools, GPU-NBDetect demonstrated several advantages. GPU-NBDetect detected the highest number of GPU-NBs among all evaluated tools, and all bugs identified by other tools were also successfully detected by GPU-NBDetect, demonstrating its comprehensive coverage. While existing tools, such as XScope and GPU-FPX, produced a significant number of false positives—17 and 33 respectively—GPU-NBDetect minimized false positives by verifying whether the generated inputs fell within the supported numerical range. Furthermore, whereas other tools often generate extreme values, such as 1.3400e+221 and 9.2579e+282, which fall outside the valid range, GPU-NBDetect systematically avoids such unsupported numbers, ensuring more reliable detection. Additionally, unlike existing approaches, GPU-NBDetect successfully detects rounding, casting, and data type bugs, which are not addressed by other tools. Verification is also more challenging with XScope, as it lacks information about the specific inputs that trigger numerical bugs.

In comparison of input generation techniques in CuPy and Numba, GPU-NBDetect outperformed RS and BO in detecting GPU-NBs and generating more effective bug-triggering inputs. Unlike RS and BO, which rely on heuristic-based sampling, GPU-NBDetect strategically mutates inputs, leading to a higher bug detection rate across Python-based GPU libraries. These findings highlight the importance of targeted input mutation techniques in uncovering hard-to-detect numerical errors. Furthermore, GPU-NBDetect achieves better overhead efficiency than existing approaches, making it more practical for large-scale GPU applications.

Some of the interesting examples observed by GPU-NBDetect are highlighted: (1) An overflow bug in cosh for input 3.6864e+4, exceeding the range specified in the CUDA Math Library, was confirmed against CPU results; (2) An underflow bug in erfc with input 6.7076092e+7 produced an output too small for floating-point representation, also verified on the CPU; and (3) A casting bug in __float2int_rd for input 1.101256458240e+11 produced the maximum 32-bit integer, indicating insufficient precision for large conversions.

Limitations of GPU-NBDetect. GPU-NBDetect has several limitations: (1) GPU-NBDetect only supports float and double data types, limiting its effectiveness for detecting errors in other numerical types. Expanding support for a broader range of data types would improve its effectiveness; (2) The tool has been evaluated on only four libraries. Future work will expand to include more GPU library APIs, but users can use GPU-NBDetect as a starting point for designing more comprehensive numerical bug detectors; (3) The mutation operators in GPU-NBDetect introduce randomness, which can lead to inconsistent detection results. Employing advanced fuzzing techniques or symbolic execution could improve the robustness of test inputs. Despite these limitations, GPU-NBDetect provides a foundational framework for test input generation and test oracle selection, paving the way for future improvements in GPU-NB detection.

Challenges in building tools to detect GPU-NBs. Developing tools to detect GPU-NBs presents unique challenges: (1) Variability in compute capabilities across GPU architectures complicates tool compatibility, limiting portability; (2) Supporting a broad range of numerical operations and data types is resource-intensive, increasing development and testing complexity; and (3) Effective GPU-NB detection requires specialized expertise in GPU computing, including compiler flags, API management, and memory allocation, which can present a steep learning curve. GPU-NBDetect provides a flexible foundation that addresses these challenges, serving as a basis for future, advanced detection tools.

Other numerical bug detection techniques. In the future, developers could improve GPU-NBDetect by adapting certain bug detection strategies for GPU-NB detection. For instance, symbolic execution [2] could be customized for GPU programs to explore diverse execution paths, helping to identify edge cases where numerical bugs may emerge. Similarly, adapting tools like Clang Static Analyzer [29] to GPU code could enable static analysis for detecting overflow, underflow, and precision errors, though this would require engineering modifications to focus on GPU-specific numerical bugs. Finally, machine learning models trained on GPU numerical bug patterns could add an advanced layer of automated detection to GPU-NBDetect, aiding in the early identification of subtle indicators of numerical instability. While each approach would need careful adaptation to suit the GPU domain, integrating these methods could significantly expand the future capabilities of GPU-NBDetect in addressing GPU-specific numerical bugs.

Fixing Numerical Bugs

Our analysis identifies common strategies developers use to manually address GPU-NBs. For each root cause, we examined common fixes based on commit histories, comments, and bug report descriptions. Table 7 summarizes these fixing strategies.

Fix	Libraries		Total
	Closed-source	Open-source
Fix through arithmetic checkers	17	10	27
Changing Precision	28	10	38
Upcast/Downcast strategies	23	9	32
Refactor Operator Implementation	39	55	94
Custom Kernel Strategy	1	4	5
Add Data Type Support/ Change Data Type	84	59	143
Fix through NaN checkers	31	16	47
Use Mixed Precision or Change Type	6	5	11

Fixing Precision Bugs

Our study outlines several strategies for addressing precision bugs in GPU programs, as shown in Figure 2. A common approach is adjusting the precision level: increasing precision can improve accuracy for out-of-range values, while decreasing precision may save memory at the cost of precision. Therefore, precision adjustments should be made with careful consideration of application requirements.

Another approach is to implement checkers that identify variables requiring precision adjustments, allowing developers to round values to appropriate levels. For example, issues like those in Figure 3 can be resolved by casting to higher precision types. Similarly, rounding bugs—such as those caused by computations with 0.5 (Figure 4)—can be mitigated by slight input adjustments (e.g., 0.49999 or 0.50001), as demonstrated in SciPy.

Casting bugs can be addressed through upcasting or downcasting. Upcasting converts values to higher precision types, while downcasting reduces precision with rounding to preserve accuracy. For example, when a Decimal64 value cannot be represented in int32 (Figure 5), casting to float64 provides a correct and effective resolution (see F-15 and I-15 in Table 1).

Fixing Data Type or Operator Bugs

Resolving data type and operator bugs in GPU programs often requires a deep understanding of CUDA numerical libraries (F-17 and I-17 in Table 1). Data type bugs are typically addressed by adding support for the required data types or substituting with compatible alternatives. Operator bugs may require refactoring numerical operations to ensure accuracy and functional correctness.

The CuPy library uses custom kernel strategies [67] such as RawKernel, ReductionKernel, and ElementwiseKernel to define operators specifically for GPU programs. These custom kernels, absent in libraries like NumPy and SciPy, are specifically designed to improve numerical stability in CuPy (F-16 and I-16 in Table 1). An example of this approach is used to fix the data type bug in Figure 9 by modifying a ReductionKernel to support operations like cupy.mean, as demonstrated in the CuPy project [62]. By explicitly defining the mean operation within ReductionKernel, users ensure more accurate computations, illustrating how specialized kernel strategies can address limitations in standard libraries.

Fixing Special Value Bugs

Fixing special value bugs involves adding checks for NaN and INF values. NaN checks identify undefined numbers, while INF checks prevent operations with values smaller than the data type’s minimum. For instance, the issue in Figure 8 can be fixed by clamping small values that are close to the numerical precision limits of float16, such as rounding values close to zero to exactly zero.

Fixing Device-specific Bugs

Fixing device-specific issues involves two main strategies: (1) switching to a GPU with better numerical support or (2) adjusting data types for mixed precision computations. For the issue in CuPy project [63], changing the data type to FP32 reduced execution time, with the NumPy version taking 0.255 seconds and the CuPy version taking 0.156 seconds, demonstrating improved performance on the GPU.

Threats to Validity

Our findings are based on a set of bugs we examined, and we applied several filtering steps to ensure correctness and relevance. First, we searched the GitHub database using strict conditions to ensure issues were related to numerical computations in GPU programming. We only selected issues confirmed as real bugs, specifically those marked as closed and containing a bug fix. We verified bug fixes by reviewing the commits associated with each bug report. Keyword filtering identified relevant numerical issues due to the large volume of issues. Each bug was independently inspected by two authors to ensure accuracy. Our results are based on both built-in and third-party numerical libraries in GPU computing, which require special attention compared to CPU computing and other GPU program bugs, enhancing generalizability.

Regarding the representativeness of our dataset, we focused on repositories associated with numerical issues in GPU programming and restricted our selection to projects impacting real applications. We employed strict criteria in our GitHub search to ensure that the samples were both targeted and representative. For methodological correctness, we collected all available commits from repositories related to numerical issues in GPU programming and utilized the GitHub Archive to find all pertinent numerical issues. We limit the numerical implications of our study to the collected dataset and focus on qualitative insights into the features of the test samples.

Related Work

CPU floating point issues vs. GPU floating point issues. Numerous studies have explored floating-point issues in CPU programming [3, 8, 11, 60, 72, 78, 80, 83, 84]. For instance, Zou et al. [84] used dynamic analysis to detect program inputs that cause floating-point errors. Wang et al. [72] proposed a method for detecting and adjusting precision in the GNU C math library. Zou et al. [83] utilized a search-based approach to generate test inputs that trigger inaccuracies. Tang et al. [60] developed a technique for obtaining high-level information about numerical stability and accuracy. Barr et al. [3] introduced Ariadne to detect floating-point exceptions in the GNU Scientific Library (GSL). Franco et al. [8] conducted a comprehensive study on numerical bugs in real-world numerical libraries.

Several studies have addressed floating-point issues in GPU programming. Laguna [26, 27, 28] introduced tools like FPCHECKER, BinFPE, and XScope to automatically detect floating-point exceptions for FP32 and FP64 data types in GPU applications. Our study focuses on numerical bugs within the supported range, addressing precision issues for various numerical operators in CUDA libraries and different data types. We incorporate differential testing to verify results against CPU outputs and use seven mutation operators for test input generation. Surineni et al. [59] examined performance-accuracy trade-offs in different arithmetic precision. Polok et al. [54] studied increasing double precision throughput on Maxwell GPUs. Ho et al. [20, 21] developed frameworks for exploiting Pascal’s half-precision capabilities and dynamically mixing precision in GPU applications. While prior studies have analyzed GPU programming bugs [17, 23, 24, 25, 31, 58, 69, 76, 79], they often focus on synchronization and memory issues rather than numerical bugs.

DL-based numerical issues vs. GPU-NBs. Prior research on numerical issues in DL has explored how flaws in algorithm design and logic, including unstable optimization methods, poorly conditioned loss functions, and incorrect gradient computations, give rise to application-level numerical instabilities [25, 53, 56, 71, 77, 78, 82]. Additionally, these studies examine the role of framework-specific transformations in exacerbating these instabilities. These issues manifest as observable deviations in model outputs, making them detectable through accuracy loss or performance degradation. In contrast, our work focuses on GPU-NBs, which arise from how data types and numerical operations are handled during execution and library implementations. These include floating-point precision limits, rounding/truncation errors, and data type mismatches. Unlike application-level issues, which impact application’s behavior, GPU-NBs silently introduce errors without explicit failure signals, affecting computations across scientific computing, high-performance computing (HPC), and DL. These issues cannot be detected through application-level techniques, necessitating low-level numerical bug detection techniques that operate independently of algorithmic design.

Detection methodologies for numerical bugs differ between DL-based approaches and GPU low-level detection techniques. DL-specific techniques, such as metamorphic and differential testing [71, 77] focus on perturbing architectures and hyperparameters to identify inconsistencies at the model level. However, these methods rely on high-level behavioral deviations, making them ineffective at exposing low-level numerical errors in GPU computations. GPU-NBs are inherently stealthy, often remaining undetected at the model’s output since their effects may not consistently propagate through multiple layers of computation.

A major limitation of DL-based approaches is the immense search space required to expose numerical bugs through model perturbation [71, 77]. Identifying arithmetic errors this way would necessitate exhaustive exploration of architectures and parameter configurations to find specific conditions where these errors amplify into observable inconsistencies. Given the combinatorial complexity of DL models, this approach is both inefficient and unreliable for detecting GPU-NBs.

In contrast, GPU-NBDetect employs a direct, low-level detection mechanism, systematically modifying floating-point precision and rounding behavior to uncover numerical inconsistencies at their source. By leveraging fine-grained bit-level perturbations and cross-device validation, GPU-NBDetect detects numerical flaws independent of model structure, addressing a fundamental gap in DL-based methodologies. The inability of DL-based testing to explicitly target numerical operations underscores the necessity of a dedicated low-level detection framework for GPU-NBs.

Device-specific numerical bugs are a common type of numerical issue, but their handling differs significantly between DL-focused studies and our approach. DL studies often attribute these bugs to library version mismatches, recommending updates as a solution [53, 56, 82]. However, this approach overlooks architectural constraints in certain GPU models, which may lack support for specific numerical operations or precision levels. Our study highlights the need for precision-aware computation strategies and effective data type selection, particularly for GPUs with limited numerical capabilities. By systematically analyzing these issues, we provide practical guidance on mitigating hardware-specific numerical errors, addressing a gap in prior research.

Beyond DL, GPU-NBs impact a wide range of GPU-based applications, including scientific computing and high-performance computing (HPC). Unlike DL numerical bugs, which primarily affect model accuracy and stability, GPU-NBs introduce fundamental computational errors that can compromise scientific results and numerical reliability in mission-critical applications. This study is the first to systematically investigate GPU-NBs at the arithmetic level. By introducing GPU-NBDetect, we provide a generalizable framework for detecting numerical inconsistencies across all GPU applications, rather than focusing solely on DL optimizations. Our findings emphasize the need for numerical bug detection techniques to ensure robust and reliable GPU-based computations.

Conclusion

This work presents a systematic investigation into GPU-NBs, offering a fine-grained characterization of their root causes, symptoms, and repair strategies. Through an empirical analysis of 397 validated bug reports from both open-source and closed-source libraries, we classify these issues into distinct categories and uncover their underlying bug patterns. Building on these insights, we propose GPU-NBDetect, a preliminary detection tool that integrates mutation-based input generation and oracle-driven validation to identify a wide range of GPU-NBs. GPU-NBDetect demonstrates strong detection capability across diverse libraries, providing a framework for automated analysis. The findings and tools introduced in this study lay the groundwork for future research on testing and repairing numerical software in GPU-based systems.

Data Availability

Our dataset and project website are available at [18].

Acknowledgments

This work was partially supported by NSF grants NSF CCF-2146443 and Amazon Research Award, Fall 2024.

References

Received 2025-02-27; accepted 2025-03-31