G.3.1 Real Vectors and Matrices
Static Semantics
1/2
2/2
generic
type Real is digits <>;
package Ada.Numerics.Generic_Real_Arrays is
pragma Pure(Generic_Real_Arrays);
3/2
-- Types
4/2
type Real_Vector is array (Integer range <>) of Real'Base;
type Real_Matrix is array (Integer range <>, Integer range <>)
of Real'Base;
5/2
-- Subprograms for Real_Vector types
6/2
-- Real_Vector arithmetic operations
7/2
function "+" (Right : Real_Vector) return Real_Vector;
function "-" (Right : Real_Vector) return Real_Vector;
function "abs" (Right : Real_Vector) return Real_Vector;
8/2
function "+" (Left, Right : Real_Vector) return Real_Vector;
function "-" (Left, Right : Real_Vector) return Real_Vector;
9/2
function "*" (Left, Right : Real_Vector) return Real'Base;
10/2
function "abs" (Right : Real_Vector) return Real'Base;
11/2
-- Real_Vector scaling operations
12/2
function "*" (Left : Real'Base; Right : Real_Vector)
return Real_Vector;
function "*" (Left : Real_Vector; Right : Real'Base)
return Real_Vector;
function "/" (Left : Real_Vector; Right : Real'Base)
return Real_Vector;
13/2
-- Other Real_Vector operations
14/2
function Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) return Real_Vector;
15/2
-- Subprograms for Real_Matrix types
16/2
-- Real_Matrix arithmetic operations
17/2
function "+" (Right : Real_Matrix) return Real_Matrix;
function "-" (Right : Real_Matrix) return Real_Matrix;
function "abs" (Right : Real_Matrix) return Real_Matrix;
function Transpose (X : Real_Matrix) return Real_Matrix;
18/2
function "+" (Left, Right : Real_Matrix) return Real_Matrix;
function "-" (Left, Right : Real_Matrix) return Real_Matrix;
function "*" (Left, Right : Real_Matrix) return Real_Matrix;
19/2
function "*" (Left, Right : Real_Vector) return Real_Matrix;
20/2
function "*" (Left : Real_Vector; Right : Real_Matrix)
return Real_Vector;
function "*" (Left : Real_Matrix; Right : Real_Vector)
return Real_Vector;
21/2
-- Real_Matrix scaling operations
22/2
function "*" (Left : Real'Base; Right : Real_Matrix)
return Real_Matrix;
function "*" (Left : Real_Matrix; Right : Real'Base)
return Real_Matrix;
function "/" (Left : Real_Matrix; Right : Real'Base)
return Real_Matrix;
23/2
-- Real_Matrix inversion and related operations
24/2
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;
25/2
-- Eigenvalues and vectors of a real symmetric matrix
26/2
function Eigenvalues (A : Real_Matrix) return Real_Vector;
27/2
procedure Eigensystem (A : in Real_Matrix;
Values : out Real_Vector;
Vectors : out Real_Matrix);
28/2
-- Other Real_Matrix operations
29/2
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1)
return Real_Matrix;
30/2
end Ada.Numerics.Generic_Real_Arrays;
31/2
{
AI95-00296-01}
The library package Numerics.Real_Arrays
is declared pure and defines the same types and subprograms as Numerics.Generic_Real_Arrays,
except that the predefined type Float is systematically substituted for
Real'Base throughout. Nongeneric equivalents for each of the other predefined
floating point types are defined similarly, with the names Numerics.Short_Real_Arrays,
Numerics.Long_Real_Arrays, etc. 31.a/2
Reason: The nongeneric
equivalents are provided to allow the programmer to construct simple
mathematical applications without being required to understand and use
generics, and to be consistent with other Numerics packages.
32/2
{
AI95-00296-01}
Two types are defined and exported by Numerics.Generic_Real_Arrays.
The composite type Real_Vector is provided to represent a vector with
components of type Real; it is defined as an unconstrained, one-dimensional
array with an index of type Integer. The composite type Real_Matrix is
provided to represent a matrix with components of type Real; it is defined
as an unconstrained, two-dimensional array with indices of type Integer.33/2
{
AI95-00296-01}
The effect of the various subprograms is as described
below. In most cases the subprograms are described in terms of corresponding
scalar operations of the type Real; any exception raised by those operations
is propagated by the array operation. Moreover, the accuracy of the result
for each individual component is as defined for the scalar operation
unless stated otherwise.34/2
{
AI95-00296-01}
In the case of those operations which are defined
to involve an inner product, Constraint_Error may be raised if
an intermediate result is outside the range of Real'Base even though
the mathematical final result would not be.35/2
function "+" (Right : Real_Vector) return Real_Vector;
function "-" (Right : Real_Vector) return Real_Vector;
function "abs" (Right : Real_Vector) return Real_Vector;
36/2
{
AI95-00296-01}
Each operation returns the result of applying the
corresponding operation of the type Real to each component of Right.
The index range of the result is Right'Range.37/2
function "+" (Left, Right : Real_Vector) return Real_Vector;
function "-" (Left, Right : Real_Vector) return Real_Vector;
38/2
{
AI95-00296-01}
Each operation returns the result of applying the
corresponding operation of the type Real to each component of Left and
the matching component of Right. The index range of the result is Left'Range.
Constraint_Error is raised if Left'Length is not equal to Right'Length.39/2
function "*" (Left, Right : Real_Vector) return Real'Base;
40/2
{
AI95-00296-01}
This operation returns the inner product of Left
and Right. Constraint_Error is raised if Left'Length is not equal to
Right'Length. This operation involves an inner product.41/2
function "abs" (Right : Real_Vector) return Real'Base;
42/2
{
AI95-00418-01}
This operation returns the L2-norm of Right (the
square root of the inner product of the vector with itself).42.a/2
Discussion: Normalization
of vectors is a frequent enough operation that it is useful to provide
the norm as a basic operation. Furthermore, implementing the norm is
not entirely straightforward, because the inner product might overflow
while the final norm does not. An implementation cannot merely return
Sqrt (X * X), it has to cope with a possible overflow of the inner product.
42.b/2
Implementation Note:
While the definition is given in terms of an inner product, the norm
doesn't “involve an inner product” in the technical sense.
The reason is that it has accuracy requirements substantially different
from those applicable to inner products; and that cancellations cannot
occur, because all the terms are positive, so there is no possibility
of intermediate overflow.
43/2
function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector;
44/2
{
AI95-00296-01}
This operation returns the result of multiplying
each component of Right by the scalar Left using the "*" operation
of the type Real. The index range of the result is Right'Range.45/2
function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector;
function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector;
46/2
{
AI95-00296-01}
Each operation returns the result of applying the
corresponding operation of the type Real to each component of Left and
to the scalar Right. The index range of the result is Left'Range.47/2
function Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) return Real_Vector;
48/2
{
AI95-00296-01}
This function returns a unit vector
with Order components and a lower bound of First. All components are
set to 0.0 except for the Index component which is set to 1.0. Constraint_Error
is raised if Index < First, Index > First + Order – 1 or
if First + Order – 1 > Integer'Last.49/2
function "+" (Right : Real_Matrix) return Real_Matrix;
function "-" (Right : Real_Matrix) return Real_Matrix;
function "abs" (Right : Real_Matrix) return Real_Matrix;
50/2
{
AI95-00296-01}
Each operation returns the result of applying the
corresponding operation of the type Real to each component of Right.
The index ranges of the result are those of Right.51/2
function Transpose (X : Real_Matrix) return Real_Matrix;
52/2
{
AI95-00296-01}
This function returns the transpose of a matrix
X. The first and second index ranges of the result are X'Range(2) and
X'Range(1) respectively.53/2
function "+" (Left, Right : Real_Matrix) return Real_Matrix;
function "-" (Left, Right : Real_Matrix) return Real_Matrix;
54/2
{
AI95-00296-01}
Each operation returns the result of applying the
corresponding operation of the type Real to each component of Left and
the matching component of Right. The index ranges of the result are those
of Left. Constraint_Error is raised if Left'Length(1) is not equal to
Right'Length(1) or Left'Length(2) is not equal to Right'Length(2).55/2
function "*" (Left, Right : Real_Matrix) return Real_Matrix;
56/2
{
AI95-00296-01}
This operation provides the standard mathematical
operation for matrix multiplication. The first and second index ranges
of the result are Left'Range(1) and Right'Range(2) respectively. Constraint_Error
is raised if Left'Length(2) is not equal to Right'Length(1). This operation
involves inner products.57/2
function "*" (Left, Right : Real_Vector) return Real_Matrix;
58/2
{
AI95-00296-01}
This operation returns the outer product of a (column)
vector Left by a (row) vector Right using the operation "*"
of the type Real for computing the individual components. The first and
second index ranges of the result are Left'Range and Right'Range respectively.59/2
function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector;
60/2
{
AI95-00296-01}
This operation provides the standard mathematical
operation for multiplication of a (row) vector Left by a matrix Right.
The index range of the (row) vector result is Right'Range(2). Constraint_Error
is raised if Left'Length is not equal to Right'Length(1). This operation
involves inner products.61/2
function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector;
62/2
{
AI95-00296-01}
This operation provides the standard mathematical
operation for multiplication of a matrix Left by a (column) vector Right.
The index range of the (column) vector result is Left'Range(1). Constraint_Error
is raised if Left'Length(2) is not equal to Right'Length. This operation
involves inner products.63/2
function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix;
64/2
{
AI95-00296-01}
This operation returns the result of multiplying
each component of Right by the scalar Left using the "*" operation
of the type Real. The index ranges of the result are those of Right.65/2
function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix;
function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix;
66/2
{
AI95-00296-01}
Each operation returns the result of applying the
corresponding operation of the type Real to each component of Left and
to the scalar Right. The index ranges of the result are those of Left.67/2
function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector;
68/2
{
AI95-00296-01}
This function returns a vector Y such that X is
(nearly) equal to A * Y. This is the standard mathematical operation
for solving a single set of linear equations. The index range of the
result is A'Range(2). Constraint_Error is raised if A'Length(1), A'Length(2),
and X'Length are not equal. Constraint_Error is raised if the matrix
A is ill-conditioned.68.a/2
Discussion: The
text says that Y is such that “X is (nearly) equal to A * Y”
rather than “X is equal to A * Y” because rounding errors
may mean that there is no value of Y such that X is exactly equal to
A * Y. On the other hand it does not mean that any old rough value will
do. The algorithm given under Implementation Advice should be followed.
68.b/2
The requirement to raise
Constraint_Error if the matrix is ill-conditioned is really a reflection
of what will happen if the matrix is ill-conditioned. See Implementation
Advice. We do not make any attempt to define ill-conditioned formally.
68.c/2
These remarks apply to
all versions of Solve and Inverse.
69/2
function Solve (A, X : Real_Matrix) return Real_Matrix;
70/2
{
AI95-00296-01}
This function returns a matrix Y such that X is
(nearly) equal to A * Y. This is the standard mathematical operation
for solving several sets of linear equations. The index ranges of the
result are A'Range(2) and X'Range(2). Constraint_Error is raised if A'Length(1),
A'Length(2), and X'Length(1) are not equal. Constraint_Error is raised
if the matrix A is ill-conditioned.71/2
function Inverse (A : Real_Matrix) return Real_Matrix;
72/2
{
AI95-00296-01}
This function returns a matrix B such that A *
B is (nearly) equal to the unit matrix. The index ranges of the result
are A'Range(2) and A'Range(1). Constraint_Error is raised if A'Length(1)
is not equal to A'Length(2). Constraint_Error is raised if the matrix
A is ill-conditioned.73/2
function Determinant (A : Real_Matrix) return Real'Base;
74/2
{
AI95-00296-01}
This function returns the determinant of the matrix
A. Constraint_Error is raised if A'Length(1) is not equal to A'Length(2).75/2
function Eigenvalues(A : Real_Matrix) return Real_Vector;
76/2
{
AI95-00296-01}
This function returns the eigenvalues of the symmetric
matrix A as a vector sorted into order with the largest first. Constraint_Error
is raised if A'Length(1) is not equal to A'Length(2). The index range
of the result is A'Range(1). Argument_Error is raised if the matrix A
is not symmetric.77/2
procedure Eigensystem(A : in Real_Matrix;
Values : out Real_Vector;
Vectors : out Real_Matrix);
78/3
{
AI95-00296-01}
{
AI05-0047-1}
This procedure computes both the eigenvalues and
eigenvectors of the symmetric matrix A. The out parameter Values is the
same as that obtained by calling the function Eigenvalues. The out parameter
Vectors is a matrix whose columns are the eigenvectors of the matrix
A. The order of the columns corresponds to the order of the eigenvalues.
The eigenvectors are normalized and mutually orthogonal (they are orthonormal),
including when there are repeated eigenvalues. Constraint_Error is raised
if A'Length(1) is not equal to A'Length(2),
or if Values'Range is not equal to A'Range(1), or if the.
The index ranges of the parameter
Vectors are not equal to those
of A. Argument_Error is raised if the matrix A is not symmetric. Constraint_Error is also raised in implementation-defined circumstances
if the algorithm used does not converge quickly enough.78.a/3
Ramification: {
AI05-0047-1}
There is no requirement on the absolute direction
of the returned eigenvectors. Thus they might be multiplied by -1. It
is only the ratios of the components that matter. This is standard practice.
79/2
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1) return Real_Matrix;
80/2
{
AI95-00296-01}
This function returns a square unit matrix
with Order**2 components and lower bounds of First_1 and First_2 (for
the first and second index ranges respectively). All components are set
to 0.0 except for the main diagonal, whose components are set to 1.0.
Constraint_Error is raised if First_1 + Order – 1 > Integer'Last
or First_2 + Order – 1 > Integer'Last.Implementation Requirements
81/2
{
AI95-00296-01}
Accuracy requirements for the subprograms Solve,
Inverse, Determinant, Eigenvalues and Eigensystem are implementation
defined. 81.a/2
Implementation defined:
The accuracy requirements for the subprograms
Solve, Inverse, Determinant, Eigenvalues and Eigensystem for type Real_Matrix.
82/2
For operations not involving
an inner product, the accuracy requirements are those of the corresponding
operations of the type Real in both the strict mode and the relaxed mode
(see G.2).
83/2
For
operations involving an inner product, no requirements are specified
in the relaxed mode. In the strict mode the modulus of the absolute error
of the inner product X*Y shall not exceed g*abs(X)*abs(Y)
where g is defined as
84/2
g = X'Length * Real'Machine_Radix**(1 – Real'Model_Mantissa)
85/2
{
AI95-00418-01}
For the L2-norm, no accuracy requirements are specified
in the relaxed mode. In the strict mode the relative error on the norm
shall not exceed g / 2.0 + 3.0 * Real'Model_Epsilon where g
is defined as above.85.a/2
Reason: This is
simply the combination of the error on the inner product with the error
on Sqrt. A first order computation would lead to 2.0 * Real'Model_Epsilon
above, but we are adding an extra Real'Model_Epsilon to account for higher
order effects.
Documentation Requirements
86/2
{
AI95-00296-01}
Implementations shall document any techniques used
to reduce cancellation errors such as extended precision arithmetic.
86.a/2
Documentation Requirement:
Any techniques used to reduce cancellation
errors in Numerics.Generic_Real_Arrays shall be documented.
86.b/2
Implementation Note:
The above accuracy requirement is met by the canonical implementation
of the inner product by multiplication and addition using the corresponding
operations of type Real'Base and performing the cumulative addition using
ascending indices. Note however, that some hardware provides special
operations for the computation of the inner product and although these
may be fast they may not meet the accuracy requirement specified. See
Accuracy and Stability of Numerical Algorithms By N J Higham (ISBN 0-89871-355-2),
Section 3.1.
86.c/3
{
AI05-0047-1}
Note moreover that the componentwise accuracy requirements
are not met by subcubic methods for matrix multiplication such as that
devised by Strassen. These methods, which are typically used for the
fast multiplication of very large matrices (e.g. order more than a few
thousands), have normwise accuracy properties. If it is desired to use
such methods, then distinct subprograms should be provided (perhaps in
a child package). See Section 22.2.2 in the above reference. Implementation Permissions
87/2
{
AI95-00296-01}
The nongeneric equivalent packages may, but need
not, be actual instantiations of the generic package for the appropriate
predefined type.Implementation Advice
88/3
{
AI95-00296-01}
{
AI05-0264-1}
Implementations should implement the Solve and
Inverse functions using established techniques such as LU decomposition
with row interchanges followed by back and forward substitution. Implementations
are recommended to refine the result by performing an iteration on the
residuals; if this is done, then it should be documented. 88.a/2
Implementation Advice:
Solve and Inverse for Numerics.Generic_Real_Arrays
should be implemented using established techniques such as LU decomposition
and the result should be refined by an iteration on the residuals.
89/2
It is not the intention that
any special provision should be made to determine whether a matrix is
ill-conditioned or not. The naturally occurring overflow (including division
by zero) which will result from executing these functions with an ill-conditioned
matrix and thus raise Constraint_Error is sufficient.
89.a/2
Discussion: There
isn't any advice for the implementation to document with this paragraph.
90/2
The test that a matrix is
symmetric should be performed by using the equality operator to compare
the relevant components.
90.a/2
Implementation Advice:
The equality operator should be used
to test that a matrix in Numerics.Generic_Real_Arrays is symmetric.
91/3
{
AI05-0047-1}
An implementation should minimize the circumstances
under which the algorithm used for Eigenvalues and Eigensystem fails
to converge. 91.a.1/3
Implementation Advice:
An implementation should minimize the
circumstances under which the algorithm used for Numerics.Generic_Real_Arrays.Eigenvalues
and Numerics.Generic_Real_Arrays.Eigensystem fails to converge.
91.a/3
Implementation Note:
J. H. Wilkinson is the acknowledged expert in this area. See for
example Wilkinson, J. H., and Reinsch, C. , Linear Algebra , vol II of
Handbook for Automatic Computation, Springer-Verlag, or Wilkinson, J.
H., The Algebraic Eigenvalue Problem, Oxford University Press.
Extensions to Ada 95
91.b/2
{
AI95-00296-01}
The package Numerics.Generic_Real_Arrays
and its nongeneric equivalents are new. Wording Changes from Ada 2005
91.c/3
{
AI05-0047-1}
Correction: Corrected various accuracy and
definition issues.
Ada 2005 and 2012 Editions sponsored in part by Ada-Europe
