G.3.1 Real Vectors and Matrices
Static Semantics
generic
type Real
is digits <>;
package Ada.Numerics.Generic_Real_Arrays
is
pragma Pure(Generic_Real_Arrays);
-- Types
type Real_Vector
is array (Integer
range <>)
of Real'Base;
type Real_Matrix
is array (Integer
range <>, Integer
range <>)
of Real'Base;
-- Subprograms for Real_Vector types
-- Real_Vector arithmetic operations
function "+" (Right : Real_Vector) return Real_Vector;
function "-" (Right : Real_Vector) return Real_Vector;
function "abs" (Right : Real_Vector) return Real_Vector;
function "+" (Left, Right : Real_Vector) return Real_Vector;
function "-" (Left, Right : Real_Vector) return Real_Vector;
function "*" (Left, Right : Real_Vector) return Real'Base;
function "abs" (Right : Real_Vector) return Real'Base;
-- Real_Vector scaling operations
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;
-- Other Real_Vector operations
function Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1)
return Real_Vector;
-- Subprograms for Real_Matrix types
-- Real_Matrix arithmetic operations
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;
function "+" (Left, Right : Real_Matrix) return Real_Matrix;
function "-" (Left, Right : Real_Matrix) return Real_Matrix;
function "*" (Left, Right : Real_Matrix) return Real_Matrix;
function "*" (Left, Right : Real_Vector) return Real_Matrix;
function "*" (Left : Real_Vector; Right : Real_Matrix)
return Real_Vector;
function "*" (Left : Real_Matrix; Right : Real_Vector)
return Real_Vector;
-- Real_Matrix scaling operations
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;
-- 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);
-- Other Real_Matrix operations
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1)
return Real_Matrix;
end Ada.Numerics.Generic_Real_Arrays;
{
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.
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.
{
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.
{
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.
{
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.
function "+" (Right : Real_Vector) return Real_Vector;
function "-" (Right : Real_Vector) return Real_Vector;
function "abs" (Right : Real_Vector) return Real_Vector;
{
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.
function "+" (Left, Right : Real_Vector) return Real_Vector;
function "-" (Left, Right : Real_Vector) return Real_Vector;
{
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.
function "*" (Left, Right : Real_Vector) return Real'Base;
{
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.
function "abs" (Right : Real_Vector) return Real'Base;
{
AI95-00418-01}
This operation returns the L2-norm of Right (the square root of the inner
product of the vector with itself).
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.
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.
function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector;
{
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.
function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector;
function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector;
{
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.
function Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) return Real_Vector;
{
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.
function "+" (Right : Real_Matrix) return Real_Matrix;
function "-" (Right : Real_Matrix) return Real_Matrix;
function "abs" (Right : Real_Matrix) return Real_Matrix;
{
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.
function Transpose (X : Real_Matrix) return Real_Matrix;
{
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.
function "+" (Left, Right : Real_Matrix) return Real_Matrix;
function "-" (Left, Right : Real_Matrix) return Real_Matrix;
{
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).
function "*" (Left, Right : Real_Matrix) return Real_Matrix;
{
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.
function "*" (Left, Right : Real_Vector) return Real_Matrix;
{
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.
function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector;
{
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.
function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector;
{
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.
function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix;
{
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.
function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix;
function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix;
{
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.
function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector;
{
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.
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.
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.
These remarks apply to all versions of Solve
and Inverse.
function Solve (A, X : Real_Matrix) return Real_Matrix;
{
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.
function Inverse (A : Real_Matrix) return Real_Matrix;
{
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.
function Determinant (A : Real_Matrix) return Real'Base;
{
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).
function Eigenvalues(A : Real_Matrix) return Real_Vector;
{
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.
procedure Eigensystem(A : in Real_Matrix;
Values : out Real_Vector;
Vectors : out Real_Matrix);
{
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 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.
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.
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1) return Real_Matrix;
{
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
{
AI95-00296-01}
Accuracy requirements for the subprograms Solve, Inverse, Determinant,
Eigenvalues and Eigensystem are implementation defined.
Implementation defined: The accuracy
requirements for the subprograms Solve, Inverse, Determinant, Eigenvalues
and Eigensystem for type Real_Matrix.
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).
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
g = X'Length * Real'Machine_Radix**(1 – Real'Model_Mantissa)
{
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.
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
{
AI95-00296-01}
Implementations shall document any techniques used to reduce cancellation
errors such as extended precision arithmetic.
Documentation Requirement: Any techniques
used to reduce cancellation errors in Numerics.Generic_Real_Arrays shall
be documented.
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.
{
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
{
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
{
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.
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.
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.
Discussion: There isn't any advice for
the implementation to document with this paragraph.
The test that a matrix is symmetric should be performed
by using the equality operator to compare the relevant components.
Implementation Advice: The equality operator
should be used to test that a matrix in Numerics.Generic_Real_Arrays
is symmetric.
{
AI05-0047-1}
An implementation should minimize the circumstances under which the algorithm
used for Eigenvalues and Eigensystem fails to converge.
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.
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
{
AI95-00296-01}
The package Numerics.Generic_Real_Arrays and its
nongeneric equivalents are new.
Wording Changes from Ada 2005
{
AI05-0047-1}
Correction: Corrected various accuracy and definition issues.
Ada 2005 and 2012 Editions sponsored in part by Ada-Europe
