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

11430

Generic package of elementary functions for Ada,

11729

Generic package of primitive functions for Ada,

13813

Generic package of real and complex type declarations and basic operations
for Ada (including vector and matrix types),

13814

Generic package of complex elementary functions for Ada.

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

a)

a complex types package including various complex arithmetic operations,

b)

a real arrays package covering both vectors and matrices,

c)

a complex arrays package covering both vectors and matrices, and

d)

a complex input–output package.

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_Vector scaling operations*

... --*operations "*" and "/" to multiply a vector by a scalar and divide a vector by a scalar*

... --

-- *Other Real_Vector operations*

**function** Unit_Vector(Index: Integer; Order: Positive; First: Integer := 1) **return** Real_Vector;

-- *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 scaling operations*

... --*operations "*" and "/" to multiply a matrix by a scalar and divide a matrix by a scalar*

... --

-- *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);

A:

Values:

-- *Other Real_Matrix operations*

**function** Unit_Matrix(Order: Positive; First_1, First_2: Integer := 1) **return** 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*_{i}*x*_{i})
– or just √(*x*_{i}*x*_{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

then we might write

X, Y: Real_Vector(1 .. N);

A: Real_Matrix(1 .. N, 1 .. N);

...

Y := Solve(A, X);

A: Real_Matrix(1 .. N, 1 .. N);

...

Y := Solve(A, X);

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

XX, YY: Real_Matrix(1 .. N, 1 .. M)

A: Real_Matrix(1 .. N, 1 .. N);

...

YY := Solve(A, XX);

A: Real_Matrix(1 .. N, 1 .. N);

...

YY := Solve(A, XX);

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

Y := Inverse(A)*X;

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^{30}.
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

to provide commonly required simple facilities
for the user who is not a numerical professional,

to provide a baseline of types and operations that
forms a firm foundation for binding to more general facilities such as
the BLAS.

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

...

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.

© 2005, 2006, 2007 John Barnes Informatics.

Sponsored in part by:

The Ada Resource Association and its member companies: |
and Ada-Europe: |