Floating-point Expansions Research Articles

Reproducibility of parallel preconditioned conjugate gradient in hybrid programming environments

The Preconditioned Conjugate Gradient method is often employed for the solution of linear systems of equations arising in numerical simulations of physical phenomena. While being widely used, the solver is also known for its lack of accuracy while computing the residual. In this article, we propose two algorithmic solutions that originate from the ExBLAS project to enhance the accuracy of the solver as well as to ensure its reproducibility in a hybrid MPI + OpenMP tasks programming environment. One is based on ExBLAS and preserves every bit of information until the final rounding, while the other relies upon floating-point expansions and, hence, expands the intermediate precision. Instead of converting the entire solver into its ExBLAS-related implementation, we identify those parts that violate reproducibility/non-associativity, secure them, and combine this with the sequential executions. These algorithmic strategies are reinforced with programmability suggestions to assure deterministic executions. Finally, we verify these approaches on two modern HPC systems: both versions deliver reproducible number of iterations, residuals, direct errors, and vector-solutions for the overhead of less than 37.7% on 768 cores.

Algorithms for Triple-Word Arithmetic

Triple-word arithmetic consists in representing high-precision numbers as the unevaluated sum of three floating-point numbers (with “nonoverlapping” constraints that are explicited in the paper). We introduce and analyze various algorithms for manipulating triple-word numbers: rounding a triple-word number to a floating-point number, adding, multiplying, dividing, and computing square-roots of triple-word numbers, etc. We compare our algorithms, implemented in the Campary library, with other solutions of comparable accuracy. It turns out that our new algorithms are significantly faster than what one would obtain by just using the usual floating-point expansion algorithms in the special case of expansions of length 3.

Numerical reproducibility for the parallel reduction on multi- and many-core architectures

A parallel algorithm to compute correctly-rounded floating-point sumsHighly-optimized implementations for modern CPUs, GPUs and Xeon PhiAs fast as memory bandwidth allows for large sums with moderate dynamic rangeScales well with the problem size and resources used on a cluster of compute nodes On modern multi-core, many-core, and heterogeneous architectures, floating-point computations, especially reductions, may become non-deterministic and, therefore, non-reproducible mainly due to the non-associativity of floating-point operations. We introduce an approach to compute the correctly rounded sums of large floating-point vectors accurately and efficiently, achieving deterministic results by construction. Our multi-level algorithm consists of two main stages: first, a filtering stage that relies on fast vectorized floating-point expansion; second, an accumulation stage based on superaccumulators in a high-radix carry-save representation. We present implementations on recent Intel desktop and server processors, Intel Xeon Phi co-processors, and both AMD and NVIDIA GPUs. We show that numerical reproducibility and bit-perfect accuracy can be achieved at no additional cost for large sums that have dynamic ranges of up to 90 orders of magnitude by leveraging arithmetic units that are left underused by standard reduction algorithms.

Arithmetic Algorithms for Extended Precision Using Floating-Point Expansions

Many numerical problems require a higher computing precision than the one offered by standard floating-point (FP) formats. One common way of extending the precision is to represent numbers in a <i>multiple component</i> format. By using the so-called <i>floating-point expansions</i> , real numbers are represented as the unevaluated sum of standard machine precision FP numbers. This representation offers the simplicity of using directly available, hardware implemented and highly optimized, FP operations. It is used by multiple-precision libraries such as Bailey's QD or the analogue Graphics Processing Units (GPU) tuned version, GQD. In this article we briefly revisit algorithms for adding and multiplying FP expansions, then we introduce and prove new algorithms for normalizing, dividing and square rooting of FP expansions. The new method used for computing the reciprocal <inline-formula><tex-math notation="LaTeX">${a}^{-1}$</tex-math></inline-formula> and the square root <inline-formula><tex-math notation="LaTeX">$\sqrt{a}$</tex-math></inline-formula> of a FP expansion <inline-formula><tex-math notation="LaTeX">$a$</tex-math></inline-formula> is based on an adapted Newton-Raphson iteration where the intermediate calculations are done using “truncated” operations (additions, multiplications) involving FP expansions. We give here a thorough error analysis showing that it allows very accurate computations. More precisely, after <inline-formula><tex-math notation="LaTeX">$q$</tex-math> </inline-formula> iterations, the computed FP expansion <inline-formula><tex-math notation="LaTeX">$x=x_0+\ldots +x_{2^q-1}$</tex-math> </inline-formula> satisfies, for the reciprocal algorithm, the relative error bound: <inline-formula><tex-math notation="LaTeX"> $\left|({x-a^{-1}})/{a^{-1}}\right| \le 2^{-2^q(p-3)-1}$</tex-math></inline-formula> and, respectively, for the square root one: <inline-formula><tex-math notation="LaTeX"> $\left|x-{1}/{\sqrt{a}}\right| \le {2^{-2^q(p-3)-1}}/{\sqrt{a}}$</tex-math></inline-formula> , where <inline-formula> <tex-math notation="LaTeX">$p&gt; 2$</tex-math></inline-formula> is the precision of the FP representation used ( <inline-formula><tex-math notation="LaTeX">$p=24$</tex-math></inline-formula> for single precision and <inline-formula><tex-math notation="LaTeX">$p=53$</tex-math></inline-formula> for double precision).