- Ada_BLAS - the non-generic parent containing definitions used by both the real and complex BLAS.
- Ada_BLAS.Real - a generic child with bindings to the single and double precision real BLAS. The precision to be used is automatically determined from the floating point type used in the instantiation. If the floating point type does not correspond to single or double precision then the routines will raise Unsupported_Precision_Error.
- Ada_BLAS.Complex - a generic child with bindings to the single and double precision complex BLAS. The precision to be used is automatically determined from the floating point type used in the instantiation. If the floating point type does not correspond to single or double precision then the routines will raise Unsupported_Precision_Error.

Please read the package specifications. You can browse the source here.

You will also find the following, which while not strictly speaking part of the binding will hopefully be helpful:

- Ada_BLAS.Get_Precision - a generic package that determines the BLAS precision corresponding to the floating point type it is instantiated with.
- Print_Precisions - a program that prints the BLAS precisions corresponding to a sample of floating point types.
- Example* - some examples.

- The matrix type used to instantiate
Ada_BLAS.Real and
Ada_BLAS.Complex
**must**follow Fortran conventions. Usepragma Convention (Fortran, Matrix_Type);

to ensure this. The same goes for the vector type used in the instantiations, though this is less critical. Usepragma Convention (Fortran, Vector_Type);

in order to sleep better at night. - If you instantiate with vectors and matrices of a constrained floating
point type, be aware that the constraints will
**not**be respected. For example, suppose you use the definitionssubtype Float_Type is Float range -1.0 .. 1.0; type Vector_Type is array (Positive range <>) of Float_Type; pragma Convention (Fortran, Vector_Type); type Matrix_Type is array (Positive range <>, Positive range <>) of Float_Type; pragma Convention (Fortran, Matrix_Type);

If you use, say, the BLAS routine AXPY to add the vectors (1.0, 1.0) and (1.0, 1.0), then it will happily return (2.0, 2.0) in spite of the fact that 2.0 is not a valid value of Float_Type. You can check the validity of returned values using the 'Valid attribute.

- Compile with optimization on in order to get dead code elimination.
All subroutines for unused floating point precisions should be eliminated
at compile time. Most subroutines should also be inlined. The only
inevitable overhead comes from the extra checking performed, and even
this can sometimes be eliminated at compile time. You should probably check
what your compiler is doing by looking at the assembly code
(use the -S switch with the GNAT compiler for example).
- LAPACK
comes with an incomplete version of the BLAS.
The missing routines can be picked up from
NETLIB.
If you are using the BLAS that come with LAPACK
and your linker is complaining that it can't find some
routines, now you know why!
- You may find that all BLAS routines for the precision
you are using are being linked, and not just the ones you
use. You may be able to overcome this by compiling your
entire program using optimization. The GNAT compiler comes
with a tool "gnatelim" that may also be helpful.
- Version 3.13 and earlier of the GNAT compiler incorrectly report an
error on line 46 of Ada_BLAS.Real. The work around is to change that line from
type Vector_Type is array (Index_Type range <>) of Float_Type'Base;

totype Vector_Type is array (Index_Type range <>) of Float_Type;

- Check that the your Ada compiler generates the right code for calling your BLAS library. For example, version 3.13 and earlier of the GNAT Ada compiler and version 3.0 and earlier of the g77 Fortran compiler differ in how to deal with arguments of type Character_Set / CHARACTER*1: the code compiled with g77 expects to be called with an extra (hidden) parameter, giving the length, for each variable of this type. The GNAT compiler does not supply this parameter, leading to occasional crashes. I am still thinking about the best solution for this problem.

SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )corresponds to

procedure GEMV ( TRANS : in Transpose_Type; M : in Natural; N : in Natural; ALPHA : in Float_Type'Base; A : in Matrix_Type; X : in Vector_Type; INCX : in Integer; BETA : in Float_Type'Base; Y : in out Vector_Type; INCY : in Integer );The translation rules are:

- The first character is dropped from the name: the Fortran subroutines SGEMV and DGEMV correspond to the single Ada subprogram GEMV. The dropped character indicates the precision of the Fortran subroutine. The binding automatically calls the Fortran subroutine of the correct precision based on the floating point type used in the instantiation (or raises Unsupported_Precision_Error). This routing to the correct subroutine should involve no cost if the Ada compiler performs dead code elimination.
- Redundant parameters are dropped. For example, DGEMV takes a parameter LDA while GEMV does not. This is because the binding can determine the correct value for LDA automatically: LDA is the leading dimension of A, which is Fortran speak for A'Length (1).
- Character parameters are replaced with enumeration types defined in
Ada_BLAS. For example,
the TRANS parameter of DGEMV is the character 'N' (do not transpose A),
'T' (transpose A) or 'C' (conjugate transpose A). The corresponding
TRANS parameter for GEMV is of Transpose_Type, defined as
type Transpose_Type is ( None, -- use X Transpose, -- use X^T Conjugate_Transpose -- use X^H );

The cost of converting between these enumeration types and Fortran Characters is small. Moreover, in the typical case that the value of the parameter is known at compile time, an optimizing compiler that performs inlining should be able to eliminate the corresponding code, resulting in no run time cost. - INTEGER arguments are converted to Integer, Natural or Positive. If an integer valued parameter must be positive, then it is declared to be of type Positive. If an integer valued parameter must be non-negative, then it is declared to be of type Natural (for example, M and N in DGEMV). Otherwise, when negative values are allowed, it is declared to be of type Integer. This seemed like a good idea at the time, but may be a mistake since Ada and Fortran integer types need not be of the same size (leading to conversion costs and possible overflow problems).
- Arguments are checked. If a check fails then
Argument_Error is raised. All constraints noted in the BLAS documentation
are checked for, as well (I hope) as all constraints which should have been
documented but weren't. For example, the DGEMV routine specifies
* LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit.

Recall that LDA is the same as A'Length (1). Thus GEMV checks that A'Length (1) is at least 1 and is not smaller than M, and raises Argument_Error if not. Checking arguments involves a small fixed cost which in some cases can be eliminated by an optimizing compiler. There are two reasons for performing these checks. The first reason is that it is in keeping with the Ada philosophy of correctness and robustness. The second reason is that the level 2 and 3 BLAS perform argument checking, and call XERBLA if a check fails. The default implementation of XERBLA prints a message and exits the program, neither of which may be desired. The only portable way I could see of replacing this with a more Ada-like mechanism was to perform the checks in the binding before calling the BLAS subroutine. This means that if a level 2 or 3 routine is called with correct parameters, then all checks are performed twice, once by the binding and once by the Fortran code. Note that the level 1 BLAS perform no error checking but that the binding does. I chose to do this in spite of the fact that the cost is more likely to be significant for level 1 routines than for level 2 or 3. - A simplified version of each routine is provided. The simplified
version assumes that all matrices and vectors are used in their entirety.
For example, the simplified version of GEMV is
procedure GEMV ( TRANS : in Transpose_Type; ALPHA : in Float_Type'Base; A : in Matrix_Type; X : in Vector_Type; BETA : in Float_Type'Base; Y : in out Vector_Type );

The missing parameters are deduced from the others: here M and N are set to A'Length (1) and A'Length (2) respectively, and INCX and INCY to 1. - Each routine is briefly documented. This is typically taken from the BLAS
quick reference guide. For GEMV it is
-- GEMV: y <- alpha A x + beta y, y <- alpha A^T x + beta y; A - M x N

I have also tried to document obscure points, for example-- DSDOT: dot <- alpha + x^T y (double precision accumulation) -- Either SDSDOT (single precision) or ALPHA + DDOT (double precision)