An Ada 95 binding to the BLAS

Introduction

The "BLAS" is a Fortran library of Basic Linear Algebra Subprograms, routines for doing fundamental vector and matrix operations. These operations are grouped into "level 1", "level 2" and "level 3". Level 1 performs vector-vector operations, for example dot products. Level 2 performs matrix-vector operations, such as the product of a matrix and a vector. Level 3 performs matrix-matrix operations, for example multiplying two matrices. Efficient implementations of the BLAS exist for many platforms. A reference implementation and documentation can be obtained from NETLIB. Thanks to their efficiency, portability and availability, the BLAS are widely used, for example by LAPACK.

Ada 95

Ada is a high level standardized (ANSI/ISO/IEC-8652:1995) programming language. Ada 95 is the latest revision of the original standard. Ada is typically used in applications where robustness and correctness are of particular importance. Information about Ada, Ada tools, compilers and the Ada community can be found at www.adapower.com.

The binding

Here is the binding as a compressed tar file. In addition to the binding, the tar file contains a copy of this page as well as the copyright notices. The binding is also available as a text file in which the various routines are simply listed one after the other. Depending on your compiler, you may need to cut it up into the individual files. A tool like gnatchop (part of the GNAT distribution) can do this for you.

Copyright

The copyright for this binding is held by the Centre National de la Recherche Scientifique (CNRS). The copyright and licence notice can be found here. Note that, as described precisely in the notice, the binding is distributed under the GNU General Public License with a special exception that allows programs to use the binding without necessarily themselves becoming subject to the GNU General Public License as a consequence.

Structure

The binding consists of the following packages:

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.

Warnings

The matrix type used to instantiate Ada_BLAS.Real and Ada_BLAS.Complex must follow Fortran conventions. Use
```
   pragma Convention (Fortran, Matrix_Type);
```
to ensure this. The same goes for the vector type used in the instantiations, though this is less critical. Use
```
   pragma 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 definitions
```
   subtype 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.

Advice and troubleshooting

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;
```
to
```
   type 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.

General principles

This is a thin binding. Each Fortran subroutine has been directly translated into an Ada subprogram. Here is an example that illustrates the rules of translation. The Fortran subroutine

      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)

Feedback

Please send me bug reports, suggestion, comments, jokes and so forth.

Back to Duncan's home page.