Table of Contents
1. Introduction
2. Content
3. Installation
4. Testing
5. Feedback
1. Introduction
This library of routines is part of a reference implementation for
the Dense and Banded BLAS routines, along with their
Extended and Mixed Precision versions, as documented in
Chapters 2 and 4 of the new BLAS Standard, which is available from:
http://www.netlib.org/blas/blast-forum/
EXTENDED PRECISION is only used internally; the input and output
arguments remain the same as in the existing BLAS. At present, we
only allow Single, Double, or Extra internal precision. Extra
precision is implemented as double-double precision (128-bit total,
106-bit significand). The routines for the double-double precision
basic arithmetic operations +, -, *, / were developed by David
Bailey.
We have designed all our routines assuming that single precision
arithmetic is actually done in IEEE single precision (32 bits)
and that double precision arithmetic is actually done in
IEEE double precision (64 bits). The routines also pass our tests
on an Intel machine with 80-bit floating point registers.
MIXED PRECISION permits some input/output arguments to be of
different types (mixing real and complex) or precisions (mixing
single and double).
The purpose of this implementation is to do a proof of
concept implementation, showing that the considerable complexity
of the specification is actually implementable and testable
with a manageable amount of software. We have not attempted to
optimize performance, but our code should be as good as
straightforward but careful code written by hand.
2. Content
The BLAS Standard defines language bindings for Fortran 95, Fortran
77, and C. Here, we have only implemented the C version and
provided a method for binding to one Fortran 77 ABI.
In this initial release, we provide the following 11 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)
SBMV (Symmetric banded matrix vector product)
SPMV (Symmetric matrix vector product, packed format)
HEMV (Hermitian matrix vector product)
HBMV (Hermitian banded matrix vector product)
HPMV (Hermitian matrix vector product, packed format)
GE_SUM_MV (Summed matrix-vector product)
TRSV (Triangular solve)
Level 3
GEMM (General matrix matrix product)
SYMM (Symmetric matrix matrix product)
HEMM (Hermitian matrix matrix product)
All have passed our systematic testing of all
possible combinations of mixed and extended precision.
We will eventually include everything in the intersection of
Chapter 2 and Chapter 4 with systematic testing.
2.1 Directory Structure
This release contains the following directories:
doc/ - technical report
src/ - Directory where C code is stored
src/dot \
src/sum } Target directories for C code
src/... /
testing/ - Directory where C code for testing is stored
testing/test-dot - DOT test code and results
testing/test-sum - SUM test code and results
testing/...
3. Installation
The reference XBLAS are built similarly to the reference BLAS and
LAPACK. You need to provide a make.inc file in the source
directory that defines the compiler, optimization flags, and
options for building a Fortran->C bridge. Some examples are
provided. The current build system produces a static libxblas.a.
The Fortran->C bridge uses details of a specific toolchain's binary
interface, in particular how the Fortran compiler mangles names.
See src/f2c-bridge.h for the available options. Most compilers can
support different name mangling schemes; be sure to use the *SAME*
naming options for all your Fortran code.
The Fortran->C bridge is included in libxblas.a. Each of the
bridge's object files matches *-f2c.o, so you can extract them with
ar if you need to share one libxblas.a between multiple Fortran
compilers. Example steps to strip the Fortran->C routines from
libxblas.a and place them in a separate libxblas-myfortran.a are as
follows:
ar t libxblas.a |fgrep -- -f2c.o | xargs ar x libxblas.a
ar ru libxblas-myfortran.a *-f2c.o
ar x libxblas.a *-f2c.o
rm *-f2c.o
4. Testing
The goal of the testing code is to validate the underlying
implementation. We must thoroughly test routines claiming to use extra
precision internally, where the test code is not allowed to declare any
extra precision variables or use any other extra precision facilities
not available to the code being tested. This requires great care in
generating test data.
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, i.e.,
| Computed_value - "True_value" | / Error_Bound
By design, the test ratio should be at most 1. A larger ratio indicates
that the computed result is either completely wrong, or not as
accurate as claimed in the specification.
The following section will discuss how we generate "good" inputs in
order to reveal the internal precisions actually used.
For details, see the paper in file doc/report.ps.
4.1 Testing DOT
DOT performs the following function:
r <- beta * r_in + alpha * (SUM_{i=1,n} x_i*y_i)
Assume that the result r_computed is computed as follows
precision eps_int internally,
precision eps_out when the final result is rounded on output
underflow threshold UN_int internally,
underflow threshold UN_out on output
and that additionally we compute a very accurate approximation r_acc with
precision eps_acc = 2^(-106) (double-double format)
underflow threshold UN_acc = 2^(-968)
Then the error bound satisfies
(*) |r_computed-r_acc| <= (n+2)(eps_int + eps_acc)*S + U + eps_out*|r_acc|
= Scale
where
S = |alpha| * (SUM_{i=1,n} |x_i|*|y_i|) + |beta|*|r_in|
U = (|alpha|*n+2)*(UN_int + UN_acc) + UN_out
Thus, we can confirm that r_computed 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_computed - r_acc| / Scale
is at most 1. Suppose that no underflow occurs, and that
eps_intX >> eps_int >= eps_acc.
where eps_intX is the internal precision actually used in some
buggy version of DOT that we want to test. Then we can expect
the test ratio to be as large as
(n+2)*(eps_intX + eps_acc)*S + eps_out*|r_acc|
test ratio ~ ----------------------------------------------
(n+2)*(eps_int + eps_acc)*S + eps_out*|r_acc|
If we can make r_acc very small, then this ratio will be roughly
test ratio = eps_intX / eps_int >> 1
which means the bug in DOT will be detected, and in fact the test ratio
will actually tell us how much internal precision we effectively used.
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.
4.2 Choosing input data and computing r_acc
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)
so that SUM_{i=1,3} x(i)*y(i) = x(2)*y(2) exactly.
4.3 Testing SPMV and GBMV
SPMV, GBMV, and many other BLAS2 routines perform the following function:
y <- beta * y + alpha * A * x
Testing it 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 the 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.
5. Feedback
Please send any comments or bug reports to extended_blas@cs.berkeley.edu.
This code was developed by
Xiaoye Li,
Jim Demmel,
David Bailey,
Yozo Hida,
Jimmy Iskandar,
Anil Kapur,
Michael Martin,
Brandon Thompson,
Teresa Tung,
Daniel Yoo
with help from Ben Wanzo, Berkat Tung, Weihua Shen, and Jason Riedy.
yozo
618b358bfd