Behavioral Non-Portability in Scientific Numeric Computing

Floating-point arithmetic is used in scientific software to perform calculations with (approximations of) real numbers. Despite successful efforts to standardize floating-point arithmetic reflected in the universally accepted IEEE 754 floating-point standard (the "Standard") , results of floating-point calculations are generally not portable across computer architectures and can in fact differ vastly.

There are a number of reasons for this phenomenon. One is the difference in floating-point hardware available on different architectures. For instance, the presence or absence of a fused multiply-add (FMA) unit significantly impacts the precision of calculating expressions of the form a b + c. Different ways of evaluating compound expressions are sanctioned by the Standard, as evaluation rules for expressions are mostly left to the programming language (unlike the result of basic arithmetic operations such as a + b).

Another reason for the differences in floating-point results is especially relevant for parallel architectures: for efficiency, complex expressions such as sums of many operands are split by the compiler into sub-expressions to be computed by individual threads; the result is in the end combined to obtain the final sum, as shown here for a sum of four arguments:

        reduction sum

Unfortunately, due to the loss of precision in floating-point arithmetic compared to real  arithmetic, many common laws of arithmetic known from high school no longer hold, in particular the associativity of addition. Different associations of summands in a long addition will therefore typically yield different results.

Most programmers are unaware of these vagaries of floating-point arithmetic. As a result, parallel scientific programs are susceptible to reliability and portability issues that can range from simple deviations in precision to changes of program control flow when moving from one architecture to another. This threat of non-portability stands in contrast to the promise made by parallel programming standards such as OpenCL for "write-once, run anywhere" functionality.

This research aims at tools and techniques to help programmers find "issues" with the use of floating-point arithmetic in parallel scientific code, specifically written using OpenCL. Specifically, the goal is to detect potential sources of reliability and portability deficiencies in such code that are due to dependencies of the floating-point behavior on the underlying (IEEE-compliant) architecture. This will have important implications for the reliability of scientific programs such as those used in biomedical imaging applications, climate modeling, and vehicle design.




We have developed a decision procedure that takes floating-point arithmetic formulas as input and decides whether they are satisfiable. Coupled with a front-end that translates verification problems over programs into floating-point decision problems, this procedure will eventually serve as the back-end of our overall tool chain. The current version of the tool is available here, with the Z3 theorem prover as an off-the-shelf back-end. This is work in progress; check the README for what works and what doesn't.

Here are some examples in Realizer's input language that demonstrate the kinds of problems the tool solves. UNKNOWN means that the back-end solver Z3 aborted with some kind of error message. TO means that the decision took more than 20mins.

Note especially the (unsatisfiable) example in red, in which the majority of rounding modes permit a solution.

Formula Real-arithm.
(x < y) /\ (y > 0.01) /\ (x = (2.0 * y - 10.0)) SAT x = -(2047.0/1048576.0)
y = (10483713.0/2097152.0)
x = -(2616197.0/262144.0)
y = (10737663.0/1073741824.0)
x = -(2616197.0/262144.0)
y = (5243.0/524288.0)
x = -(2616197.0/262144.0)
y = (10737419.0/1073741824.0)
Above: simple linear constraints, containing a constant that needs to be rounded. Satisfiable in the Reals and in single precision FP arithmetic
(y = x * x) /\ (x = (y + 0.5)) UNSAT UNSAT UNSAT UNSAT UNSAT UNSAT
Above: a quadratic equation in x that has no solution in the Reals or in FP numbers (roots are complex numbers).
(3.1 * x) > (y + 22.34) SAT x = (10824011.0/1048576.0)
y = -(16777215.0/1048576.0)
Above: a linear constraint where co-efficient of one variable is a constant that requires rounding. The other constant also needs rounding.
(x + (y + z)) > ((x + y) + z) UNSAT x = (8589935.0/8589934592.0)
y = -(1061159.0/268435456.0)
z = (8388607.0/4294967296.0)
x = -(8388607.0/2147483648.0)
y = -(8589939.0/8589934592.0)
z = (12683577.0/4294967296.0)
x = (8388609.0/2147483648.0)
y = -(12482247.0/8589934592.0)
z = -(12482249.0/8589934592.0)
Above: an inequality that is unsatisfiable in the Reals but has a satisfying assignment in FP arithmetic (FP addition is not associative).
(x < 3.0) /\ (((x * x) * x) >=39.28) /\ (x > -100.0) UNSAT UNSAT UNSAT UNSAT UNSAT UNSAT
Above: constraints that are unsatisfiable in the Reals. Also contains a cubic expression.
(6.875 * x) < (y + 39.0125) SAT x = - (19065.0/16777216.0))
y = (8589935.0/8589934592.0)
x = (8589935.0/8589934592.0)
y = - (511.0/262144.0))
x = - (16777215.0/8589934592.0)
y = - (511.0/262144.0)
UNKNOWN x = (9761289.0/4294967296.0)
y = - (33.0/32768.0)
Above: a simple expression involving multiplication by a constant. Both constants in the constraints can be exactly represented in single precision and hence don't need to be rounded
((x + 1.97) > 11.73) /\ (x = (y - 23.56)) SAT x = (7355761.0/524288.0)
y = (9853993.0/262144.0)
x = (15233187.0/65536.0)
y = (16777215.0/65536.0)
x = (5300551.0/131072.0)
y = (16777215.0/262144.0)
x = (8388607.0/524288.0)
y = (648151.0/16384.0)
x = (8130397.0/131072.0)
y = (11218453.0/131072.0)
Above: a conjunction of two simple linear constraints. All 3 constants that appear need to be rounded.

The tool was developed mostly by Jaideep Ramachandran.

Abstraction and Refinement based Model Lifting for Realizer

To make Realizer scalable, we are currently working on a sound procedure that much more aggressively exploits the proximity of Floating-Point Arithmetic (FPA) to the proxy theory of real arithmetic, by lifting an encountered Real Arithmetic (RA) model to a numerically close FPA model. Should this fail, our procedure lazily reduces the gap between FPA and RA interpretations of the input formula f, by encoding according to FPA semantics those parts of f that are likely responsible for the earlier spurious model. This process is iterated until an FPA model is found, or the entire formula is encoded according to FPA semantics, at which time we resort to existing bit-precise but expensive FPA solvers.



    National Science Foundation, under award number CCF-1218075.

        Science Foundation

External links (follow at your own risk!):

Thomas Wahl, 2013