A Reference Implementation for Extended and Mixed Precision BLAS

In this release, we provide the following routines:

• Level 1
  • DOT (inner product)
  • SUM (sum)
  • AXPBY (Scaled vector accumulation)
  • WAXPBY (Scaled vector addition)
• Level 2
  • GEMV (General matrix vector product)
  • GBMV (Banded matrix vector product)
  • SYMV (Symmetric matrix vector product)
  • SPMV (Symmetric matrix vector product, packed format)
  • HEMV (Hermitian matrix vector product)
  • HPMV (Hermitian matrix vector product, packed format)
  • TRSV (Triangular solve)
• Level 3
  • GEMM (General matrix matrix product)
  • SYMM (Symmetric matrix matrix product)
  • HEMM (Hermitian matrix matrix product)

For each BLAS routine, we perform the following steps in the test code:

1. Generate input scalars, vectors and arrays, according to the routine's specification, so that the result exposes the internal precision actually used.
2. Call the BLAS routine
3. For each output, compute a test ratio of the computed error to the theoretical error bound:

   	test ratio = | Computed_value - "True_value" | / Error_Bound
   

The following section will discuss how we generate "good" inputs in order to reveal the internal precisions actually used.

	r <- alpha * SUM_{i=1,n}(x_i*y_i) + beta * r_in

(*) |r_comp - r_acc| <= (n+2)*(eps_int + eps_acc)*S + U + eps_out*|r_acc|
                      = Scale

    S = |alpha| * SUM_{i=1,n} |x_i*y_i| + |beta * r_in|
    U = (2*|alpha|*n + 3)*(UN_int + UN_acc) + UN_out	

Thus, we can confirm that r_comp has been computed as accurately as claimed (i.e. with internal precision defined by eps_int and UN_int) by testing whether the

          test ratio = |r_comp - r_acc| / Scale

          eps_int_actual >> eps_int >= eps_acc.

                        (n+2)*(eps_int_actual + eps_acc)*S + eps_out*|r_acc|
          test ratio ~  ----------------------------------------------
                        (n+2)*(eps_int  + eps_acc)*S + eps_out*|r_acc|

          test ratio =  eps_int_actual / eps_int >> 1

Thus our goal is to pick test data alpha, x(1:n), y(1:n), beta and r_in to make |r_acc| as small as possible, MUCH smaller than S, so that eps_int term dominates on the right of the inequality (*). Otherwise eps_int will be invisible in the error bound, then we cannot tell what precision is used internally.

In our test generator, we choose input data alpha, beta, r_in, x(1:n) and y(1:n) judiciously in order to cause as much cancellation in r as possible.

The general approach is to choose some of the input values of x(i) and y(i) so that the exact (partial) dot product of these values has a lot of nonzero fraction bits, preferably at least 106. Then the remaining values of x(i) and y(i) are chosen to cancel the previous bits as much as possible. This latter computation seems to require high precision.

One possibility to use an arbitrary precision package, such as MPFUN, but our goal is to make this test code self contained, and use no extra precision facilities not available to the code being tested. Since most of our routines can be reduced to a series of dot products, testing can be based on DOT. We only need a TRUSTED dot routine. Currently, we are using our own dot routine with double-double internal precision to compute r_truth, which is accurate to 106 bits. This means that any internal precision higher than double-double cannot be detected, and may result in a tiny test ratio. A very tiny test ratio (such as zero) may also occur if the result happens to be computed exactly.

This raises the possibility that there are "matching" bugs in our trusted DOT and the DOT under test, so that bugs are missed. To avoid this possibility we also generate some test examples where the cancellation is done mathematically (and so exactly) rather than depending on a computation. The idea is simple: For example, choose x(1:3) and y(1:3) so that

            x(1)*y(1) = -x(3)*y(3) >> x(2)*y(2)

	y <- alpha * A * x + beta * y

Testing them is no more difficult than testing DOT, because each component of the computed y vector is a dot product, and satisfies the error bound (*). So we simply use the same test generator as DOT, and compute a test ratio for each component of y vector.The only tricky part is that some entries in each dot product may be fixed. For example, the first row of A and the vector x can be chosen freely, but after that x is fixed and, if A is symmetric, the first entry of each subsequent row is fixed. Our dot-product test generator handles all these cases.

This approach can be generalized to most other Level 2 and 3 BLAS, except for the triangular solve function, which is future work.

Note that conventional operator overloading where the implementation of each binary operation only depends on the argument types and not the result type is not enough for our purposes. It cannot implement extra = double*double correctly, for example. The alternative would be to promote one argument to the output type, but could be inefficient depending on how it is implemented (extra = extra*double is significantly more expensive than extra = double*double). This may limit the utility of this feature in F95 and C++.

Second, we will evaluate the overall performance of the library. Our code in the present form is meant to be a reference implementation, serving much the same purpose as the existing Fortran BLAS distribution in Netlib. Our generated code consists only of straightforward loops, without any architectural optimizations, such as blocking for locality or loop unrolling. Because of the high degree of data reuse in double-double arithmetic, this yielded good performance for extended precision, but further improvements are possible, especially for mixed precision.

And third, we will investigate the benefit of the new BLAS for more linear algebra algorithms and applications.