Numerics annex

When Ada 95 was being designed, the Numerics Rapporteur Group pontificated at length over what features should be included in Ada 95 itself, what should be placed in secondary standards, and what should be left to the creativeness of the user community.

A number of secondary standards had been developed for Ada 83. They were

The first two, 11430 and 11729, were incorporated into the Ada 95 core language. The elementary functions, 11430, (Sqrt, Sin, Cos etc) became the package Ada.Numerics.Generic_Elementary_ Functions in A.5.1, and the primitive functions, 11729, became the various attributes such as Floor, Ceiling, Exponent and Fraction in A.5.3. The original standards were withdrawn long ago.

The other two standards, although originally developed as Ada 83 standards did not become finally approved until 1998.

In the case of 13814, the functionality was all incorporated into the Numerics annex of Ada 95 as the package Ada.Numerics.Generic_Complex_Elementary_Functions in G.1.2. Accordingly the original standard has now lapsed.

However, the situation regarding 13813 was not so clear. It covered four areas

The first of these was incorporated into the Numerics annex of Ada 95 as the package Ada.Numerics.Generic_Complex_Types in G.1.1, and the last similarly became the package Ada.Text_IO.Complex_IO in G.1.3. However, the array packages, both real and complex, were not incorporated into Ada 95.

The reason for this omission is explained in Section G.1.1 of the Rationale for Ada 95 [3] which says

A decision was made to abbreviate the Ada 95 packages by omitting the vector and matrix types and operations. One reason was that such types and operations were largely self-evident, so that little real help would be provided by defining them in the language. Another reason was that a future version of Ada might add enhancements for array manipulation and so it would be inappropriate to lock in such operations permanently.

The sort of enhancements that perhaps were being anticipated were facilities for manipulating arbitrary subpartitions of arrays such as were provided in Algol 68. These rather specialized facilities have not been added to Ada 2005 and indeed it seems most unlikely that they would ever be added. The second justification for omitting the vector and matrix facilities of 13813 thus disappears.

In order to overcome the objection that everything is self-evident we have taken the approach that we should further add some basic facilities that are commonly required, not completely trivial to implement, but on the other hand are mathematically well understood.

So the outcome is that Ada 2005 includes almost everything from 13813 plus subprograms for

A small number of operations that were not related to linear algebra were removed (such as raising all elements of a matrix to a given power).

So Ada 2005 includes two new packages which are Ada.Numerics.Generic_Real_Arrays and Ada.Numerics.Generic_Complex_Arrays. It would take too much space to give the specifications of both in full so we give just an abbreviated form of the real package in which the specifications of the usual operators are omitted thus

   -- Types
   type Real_Vector is array (Integer range <>) of Real'Base;
   type Real_Matrix is array (Integer range <>, Integer range <>) of Real'Base;

   -- Real_Vector arithmetic operations
   ... -- unary and binary "+" and "–" giving a vector
   ... -- also inner product and two versions of "abs" – one returns a vector and the
   ... -- other a value of Real'Base

   -- Real_Matrix arithmetic operations
   ... -- unary "+", "–", "abs", binary "+", "–" giving a matrix
   ... -- "*" on two matrices giving a matrix, on a vector and a matrix giving a vector,
   ... -- outer product of two vectors giving a matrix, and of course
   function Transpose(X: Real_Matrix) return Real_Matrix;

   -- Real_Matrix inversion and related operations
   function Solve(A: Real_Matrix; X: Real_Vector) return Real_Vector;
   function Solve(A, X: Real_Matrix) return Real_Matrix;
   function Inverse(A: Real_Matrix) return Real_Matrix;
   function Determinant(A: Real_Matrix) return Real'Base;

   -- Eigenvalues and vectors of a real symmetric matrix
   function Eigenvalues(A: Real_Matrix) return Real_Vector;
   procedure Eigensystem(
         A: in Real_Matrix;
         Values: out Real_Vector; Vectors: out Real_Matrix);

Many of these operations are quite self-evident. The general idea as far as the usual arithmetic operations are concerned is that we just write an expression in the normal way as illustrated in the Introduction. But the following points should be noted.

There are two operations "abs" applying to a Real_Vector thus

One returns a vector each of whose elements is the absolute value of the corresponding element of the parameter (rather boring) and the other returns a scalar which is the so-called L2-norm of the vector. This is the square root of the inner product of the vector with itself or √(Σx_ix_i) – or just √(x_ix_i) using the summation convention (which will be familiar to those who dabble in the relative world of tensors). This is provided as a distinct operation in order to avoid any intermediate overflow that might occur if the user were to compute it directly using the inner product "*".

There are two functions Solve for solving one and several sets of linear equations respectively. Thus if we have the single set of n equations

and if we have m sets of n equations we might write

The functions Inverse and Determinant are provided for completeness although they should be used with care. Remember that it is foolish to solve a set of equations by writing

because it is both slow and prone to errors. The main problem with Determinant is that it is liable to overflow or underflow even for moderate sized matrices. Thus if the elements are of the order of a thousand and the matrix has order 10, then the magnitude of the determinant will be of the order of 10³⁰. The user may therefore have to scale the data.

Two subprograms are provided for determining the eigenvalues and eigenvectors of a symmetric matrix. These are commonly required in many calculations in domains such as elasticity, moments of inertia, confidence regions and so on. The function Eigenvalues returns the eigenvalues (which will be non-negative) as a vector with them in decreasing order. The procedure Eigensystem computes both eigenvalues and vectors; the parameter Values is the same as that obtained by calling the function Eigenvalues and the parameter Vectors is a matrix whose columns are the corresponding eigenvectors in the same order. The eigenvectors are mutually orthonormal (that is, of unit length and mutually orthogonal) even when there are repeated eigenvalues. These subprograms apply only to symmetric matrices and if the matrix is not symmetric then Argument_Error is raised.

Other errors such as the mismatch of array bounds raise Constraint_Error by analogy with built-in array operations.

The reader will observe that the facilities provided here are rather humble and presented in a simple black-box style. It is important to appreciate that we do not see the Ada predefined numerics library as being in any way in competition with or as a substitute for professional libraries such as the renowned BLAS (Basic Linear Algebra Subprograms, see www.netlib.org/blas). Indeed our overall goal is twofold

We do not expect users to apply the operations in our Ada packages to the huge matrices that arise in areas such as partial differential equations. Such matrices are often of a special nature such as banded and need the facilities of a comprehensive numerical library. We have instead striven to provide easy to use facilities for the programmer who has a small number of equations to solve such as might arise in navigational applications.

Simplicity is evident in that functions such as Solve do not reveal the almost inevitable underlying LU decomposition or provide parameters controlling for example whether additional iterations should be applied. However, implementations are advised to apply an additional iteration and should document whether they do or not.

Considerations of simplicity also led to the decision not to provide automatic scaling for the determinant or to provide functions for just the largest eigenvalue and so on.

Similarly we only provide for the eigensystems of symmetric real matrices. These are the ones that commonly arise and are well behaved. General nonsymmetric matrices can be troublesome.

Appropriate accuracy requirements are specified for the inner product and L2-norm operations. Accuracy requirements for Solve, Inverse, Determinant, Eigenvalues and Eigenvectors are implementation defined which means that the implementation must document them.

The complex package is very similar and will not be described in detail. However, the generic formal parameters are interesting. They are

with Ada.Numerics.Generic_Real_Arrays, Ada.Numerics.Generic_Complex_Types;
generic
   with package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays(<>);
   use Real_Arrays;
   with package Complex_Types is new Ada.Numerics.Generic_Complex_Types(Real);
   use Complex_Types;
package Ada.Numerics.Generic_Complex_Arrays is
   ...

Thus we see that it has two formal packages which are the corresponding real array package and the existing Ada 95 complex types and operations package. The formal parameter of the first is <> and that of the second is Real which is exported from the first package and ensures that both are instantiated with the same floating point type.

As well as the obvious array and matrix operations, the complex package also has operations for composing complex arrays from cartesian and polar real arrays, and computing the conjugate array by analogy with scalar operations in the complex types package. There are also mixed real and complex array operations but not mixed imaginary, real and complex array operations. Altogether the complex array package declares some 80 subprograms (there are around 30 in the real array package) and adding imaginary array operations would have made the package unwieldy (and the reference manual too heavy).

By analogy with real symmetric matrices, the complex package has subprograms for determining the eigensystems of Hermitian matrices. A Hermitian matrix is one whose complex conjugate equals its transpose; such matrices have real eigenvalues and are well behaved.

We conclude this discussion of the Numerics annex by mentioning one minute change regarding complex input–output. Ada 2005 includes preinstantiated forms of Ada.Text_IO.Complex_IO such as Ada.Complex_Text_IO (for when the underlying real type is the type Float), Ada.Long_Complex_Text_IO (for type Long_Float) and so on. These are by analogy with Float_Text_IO, Long_Float_Text_IO and their omission from Ada 95 was probably an oversight.

Rationale for Ada 2005

7.6 Numerics annex