GitHub

API - NumericalTest

The rigour of the Float32 and Float64 paths rests on two properties of the machine the package is running on: that every BLAS thread takes the rounding mode of its caller, and that the loops Julia emits for a sum, a difference or a scalar multiple of arrays round as the mode says. Neither is guaranteed, so both are tested when the package is loaded, and the outcome of each is recorded in a flag that the product and sum routines read. Where a test failed the operations concerned are formed entry by entry with the directed operations instead, or the product is enclosed without a rounding mode at all.

The rounding-mode tests

BallArithmetic.NumericalTest.rounding_test — Function
rounding_test(n, k)

Let u=fill(2^(-53), k-1) and let A be the matrix [I u; 0 2^(-53)]

This test checks the result of A*A' in different rounding modes, running BLAS on n threads: the diagonal entries but the last are 1 + u² exactly, so they must be nextfloat(1.0) under RoundUp and 1.0 under RoundDown.

Returns true only when both checks hold, so a false is a failure of the BLAS to carry the rounding mode across its threads and not a statement about the matrix.

This is the example of Section 5.1.3 of the thesis below (Listing 5.2), which reports that OpenBLAS 0.1.0 with more than one thread returned 1 under RoundUp in the last three quarters of the diagonal for n ≥ 16.

Reference

[36], doctoral thesis, École normale supérieure de Lyon, NNT 2014ENSL0941, https://theses.hal.science/tel-01126973.

source
BallArithmetic.NumericalTest.array_rounding_test — Function
array_rounding_test(T = Float64; n = 48) -> Vector{String}

The array operations that the package performs under a rounding mode, each compared entry by entry with the same operation done with the directed operations of RoundingEmulator. Returns the names of the forms that differ; an empty vector when all agree.

A product of floating-point matrices is an OpenBLAS call, and the build loaded by the package makes every BLAS thread take the rounding mode of the caller (rounding_test checks that). A sum of matrices, a multiple, a difference are Julia's own loops, and the OpenBLAS build has no part in them: that they round as the mode says is a property of the code the Julia compiler emits for this machine, which Julia does not document as guaranteed. This test is that check, and it is run when the package is loaded.

The forms: the radius of a sum, ϵ|C| + r_A + r_B upward; P + R upward and P − R downward, with matrix operations and with broadcasting; the centre (L + U)/2 and the radius M − L of two bounds, upward; each on dense matrices, on Symmetric ones, on Diagonal ones, and with C and P dense and the other two Diagonal (the sums that a product with a Diagonal factor leads to); and the products DC and CD with D diagonal, upward and downward.

source

What the tests recorded

BallArithmetic.NumericalTest.BLAS_ROUNDING_OK — Constant
BLAS_ROUNDING_OK

Whether rounding_test passed when the package was loaded, on all the BLAS threads or, failing that, on one. When it is false the products of Float32 and Float64 matrices that the package bounds from below and from above are not evaluated under a rounding mode: they are enclosed by the splitting scheme of [16] (INTLAB's prodK), with the Dot2 accumulation of [11] where the scheme declines, which do not use the rounding mode.

source