!standard G (00) 03-05-22 AI95-00296/02

!class amendment 02-06-07

!status work item 03-01-23

!status received 02-06-07

!priority Medium

!difficulty Medium

!subject Vector and matrix operations

!class amendment 02-06-07

!status work item 03-01-23

!status received 02-06-07

!priority Medium

!difficulty Medium

!subject Vector and matrix operations

!summary

The vector and matrix operations in ISO/IEC 13813 plus related operations
are added to Ada.Numerics in Annex G.

!problem

A number of secondary standards for Ada 83 were produced for the numerics area.
Most of the functionality of these standards was incorporated into Ada 95 (some
in the core language and some in the Numerics Annex) but two packages from
ISO/IEC 13813 were not. These are generic packages for the manipulation of real
and complex vectors and matrices.

The UK was asked to review the situation and to make a recommendation; they
recommended that if Ada were amended then consideration should be given to
including the packages within the Numerics annex.

The packages can be implemented entirely in Ada and thus present little burden
to implementors. Providing secondary standards has not proved satisfactory
because they are not sufficiently visible to the user community as a whole.

!proposal

It is proposed that two generic packages be added to the numerics annex. They
are Ada.Numerics.Generic_Real_Arrays and Ada.Numerics.Generic_Complex_Arrays.
They are included as a new subclause G.3 in order to avoid excessive
renumbering of other clauses.

In addition to the main facilities of 13813, these packages also include
operations for matrix inversion and determinants. Furthermore, child packages
provide subprograms for the determination of eigenvalues and eigenvectors.

!wording

Add a new clause G.3 as follows

G.3 Vector and matrix manipulation

Types and operations for the manipulation of real vectors and matrices are
provided in Generic_Real_Arrays, which is defined in G.3.1. Types and
operations for the manipulation of complex vectors and matrices are provided in
Generic_Complex_Arrays, which is defined in G.3.2. Both of these library units
are generic children of the predefined package Numerics (see A.5). Nongeneric
equivalents of these packages for each of the predefined floating point types
are also provided as children of Numerics.

G.3.1 Real Vectors and Matrices

The generic library package Numerics.Generic_Real_Arrays has the following
declaration:

-- Types

-- Real_Vector additive operations

-- Real_Vector scaling operations

-- Inner, Outer and Vector products

-- Other Real_Vector operations

-- Real_Matrix arithmetic operations

-- Real_Matrix scaling operations

-- Real_Matrix inversion and related operations

-- Other Real_Matrix operations

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.

Two types are defined and exported by Ada.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.

The effect of the various functions is as described below. In most cases the
functions 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.

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;

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;

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. The exception Constraint_Error
is raised if Left'Length is not equal to Right'Length.

function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector;

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;

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 "*" (Left, Right : Real_Vector) return Real'Base;

This operation returns the inner product of Left and Right. The exception
Constraint_Error is raised if Left'Length is not equal to Right'Length. This
operation involves an inner product.

function "*" (Left, Right : Real_Vector) return Real_Matrix;

This operation returns the outer product of a (column) vector Left by a
(row) vector Right using the appropriate operation "*" of the type Real for
computing the individual components. The first and second index ranges of the
matrix result are Left'Range and Right'Range respectively.

function "*" (Left, Right : Real_Vector) return Real_Vector;

This operation returns the vector product of Left and Right. The exception
Constraint_Error is raised if Left'Length and Right'Length are not equal to 3.
The index range of the result is Left'Range.

?? I am not sure about vector product. It only exists in three dimensions and
is really a skew-symmetric tensor.

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. The exception 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;

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.

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;

Each operation returns the result of applying the corresponding operation of
type Real to each component of Left and the matching component of Right.
The index ranges of the result are those of Left. The exception
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;

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. The exception
Constraint_Error is raised if Left'Length(2) is not equal to Right'Length(1).
This operation involves inner products.

function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector;

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). The exception 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;

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). The exception 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;

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
matrix 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;

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 matrix result are those of Left.

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 X'Range. 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.

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 those of X. 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.

This function returns a matrix B such that A * B is (nearly) the unit matrix.
The index ranges of the result are those of A. 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.

This function returns the determinant of the matrix A. Constraint_Error is
raised if A'Length(1) is not equal to A'Length(2).

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. The exception Constraint_Error is raised if
First_1 + Order - 1 > Integer'Last or First_2 + Order - 1 > Integer'Last.

Implementation Requirements

Accuracy requirements for the functions Solve, Inverse and Determinant are
implementation defined.

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).

?? This doesn't really cover vector product but I might delete it anyway.

For operations involving an inner product, no requirements are specified in
the relaxed mode. In the strict mode the acuracy should be at least that of
the canonical implementation of multiplication and addition using the
corresponding operations of type Real'Base and performing the cumulative
addition using ascending indices. Implementations shall document any
techniques used to reduce cancellation errors such as extended precision
arithmetic.

Implementation Permissions

The nongeneric equivalent packages may, but need not, be actual
instantiations of the generic package for the appropriate predefined type.

Implementation Advice

Implementations should implement the Solve and Inverse functions using
established techniques such as Gauss-Seidel with Interchanges.

It is not the intention that any special provision should be made to
determine whether a matrix is ill-conditioned or not. The naturally occuring
division by zero or overflow which will result from executing these
functions with an ill-conditioned matrix and thus raise Constraint_Error is
sufficient.

G.3.2 Complex Vectors and Matrices

The generic library package Numerics.Generic_Complex_Arrays has the following
declaration:

-- Types

-- Complex_Vector selection, conversion and composition operations

-- Complex_Vector additive operations

-- Mixed Real_Vector and Complex_Vector additive operations

-- Complex_Vector scaling operations

-- Inner, Outer and Vector products

-- Other Complex_Vector operations

-- Complex_Matrix selection, conversion and composition operations

-- Complex_Matrix arithmetic operations

-- Mixed Real_Matrix and Complex_Matrix arithmetic operations

-- Complex_Matrix scaling operations

-- Complex_Matrix inversion and related operations

?? I have not added mixed real and complex for Solve.

-- Other Complex_Matrix operations

function Unit_Matrix (Order : Positive;

First_1, First_2 : Integer := 1)

return Complex_Matrix;

The library package Numerics.Complex_Arrays is declared pure and defines
the same types and subprograms as Numerics.Generic_Complex_Arrays, except
that the predefined type Float is systematically substituted for Real'Base,
and the Real_Vector and Real_Matrix types exported by Numerics.Real_Arrays
are systematically substituted for Real_Vector and Real_Matrix, and the
Complex and Imaginary types exported by Numerics.Complex_Types are
systematically substituted for Complex and Imaginary, throughout. Nongeneric
equivalents for each of the other predefined floating point types are
defined similarly, with the names Numerics.Short_Complex_Arrays,
Numerics.Long_Complex_Arrays, etc.

?? I based this text on that in G.1.2(9/1). However this clause (ie G.3.2)
doesn't seem to use Imaginary but I left it in.

Two types are defined and exported by Ada.Numerics.Generic_Complex_Arrays.
The composite type Complex_Vector is provided to represent a vector with
components of type Complex; it is defined as an unconstrained
one-dimensional array with an index of type Integer. The composite type
Complex_Matrix is provided to represent a matrix with components of type
Complex; it is defined as an unconstrained, two-dimensional array with
indices of type Integer.

The effect of the various subprograms is as described below. In many cases
they are described in terms of corresponding scalar operations in
Numerics.Generic_Complex_Types. Any exception raised by those operations is
propagated by the array subprogram. Moreover any constraints on the parameters
and the accuracy of the result for each individual component is as defined for
the scalar operation.

In the case of those operations which are defined to involve an inner product,
Constraint_Error may be raised if an intermediate result has a component
outside the range of Real'Base even though the final mathematical result
would not.

Each function returns a vector of the specified cartesian components of X.
The index range of the result is X'Range.

Each procedure replaces the specified (cartesian) component of each of
the components of X by the value of the matching component of Re or Im;
the other (cartesian) component of each of the components is unchanged.
The exception Constraint_Error is raised if X'Length is not equal to
Re'Length and if X'Length is not equal to Im'Length.

Each function constructs a vector of Complex results (in cartesian
representation) formed from given vectors of cartesian components; when
only the real components are given, imaginary components of zero are assumed.
The index range of the result is Re'Range. The exception Constraint_Error is
raised if Re'Length is not equal to Im'Length.

Each function calculates and returns a vector of the specified polar
components of X or Right using the corresponding function in
Numerics.Generic_Complex_Types. The index range of the result is X'Range or
Right'Range.

Each function constructs a vector of Complex results (in cartesian
representation) formed from given vectors of polar components using the
corresponding function in Numerics.Generic_Complex_Types on matching
components of Modulus and Argument. The index range of the result is
Modulus'Range. The exception Constraint_Error is raised if
Modulus'Length is not equal to Argument'Length.

function "+" (Right : Complex_Vector) return Complex_Vector;
function "-" (Right : Complex_Vector) return Complex_Vector;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of Right. The index
range of the result is Right'Range.

This function returns the result of applying the appropriate function
Conjugate in Numerics.Generic_Complex_Types to each component of X. The
index range of the result is X'Range.

function "+" (Left, Right : Complex_Vector) return Complex_Vector;
function "-" (Left, Right : Complex_Vector) return Complex_Vector;

function "+" (Left : Real_Vector;

Right : Complex_Vector) **return** Complex_Vector;

function "+" (Left : Complex_Vector;

Right : Real_Vector) **return** Complex_Vector;

function "-" (Left : Real_Vector;

Right : Complex_Vector) **return** Complex_Vector;

function "-" (Left : Complex_Vector;

Right : Real_Vector) **return** Complex_Vector;

Each operation returns the result of applying the corresponding operation
in Numerics.Generic_Complex_Types to each component of Left and the
matching component of Right. The index range of the result is Left'Range.
The exception Constraint_Error is raised if Left'Length is not equal to
Right'Length.

function "*" (Left : Complex; Right : Complex_Vector) return Complex_Vector;

This operation returns the result of multiplying each component of Right by
the complex number Left using the appropriate operation "*" in
Numerics.Generic_Complex_Types. The index range of the result is Right'Range.

function "*" (Left : Complex_Vector; Right : Complex) return Complex_Vector;
function "/" (Left : Complex_Vector; Right : Complex) return Complex_Vector;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of the vector Left and the
complex number Right. The index range of the result is Left'Range.

function "*" (Left : Real'Base; Right : Complex_Vector) return Complex_Vector;

This operation returns the result of multiplying each component of Right by
the real number Left using the appropriate operation "*" in
Numerics.Generic_Complex_Types. The index range of the result is Right'Range.

function "*" (Left : Complex_Vector; Right : Real'Base) return Complex_Vector;
function "/" (Left : Complex_Vector; Right : Real'Base) return Complex_Vector;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of the vector Left and the
real number Right. The index range of the result is Left'Range.

function "*" (Left, Right : Complex_Vector) return Complex;
function "*" (Left : Real_Vector; Right : Complex_Vector) return Complex;
function "*" (Left : Complex_Vector; Right : Real_Vector) return Complex;

Each operation returns the inner product of Left and Right. The exception
Constraint_Error is raised if Left'Length is not equal to Right'Length.
These operations involve an inner product.

function "*" (Left, Right : Complex_Vector) return Complex_Matrix;
function "*" (Left : Real_Vector;

Right : Complex_Vector) **return** Complex_Matrix;

function "*" (Left : Complex_Vector;

Right : Real_Vector) **return** Complex_Matrix;

Each operation returns the outer product of a (column) vector Left by a (row)
vector Right using the appropriate operation "*" in
Numerics.Generic_Complex_Types for computing the individual components.
The first and second index ranges of the matrix result are Left'Range and
Right'Range respectively.

function "*" (Left, Right : Complex_Vector) return Complex_Vector;
function "*" (Left : Real_Vector;

Right : Complex_Vector) **return** Complex_Vector;

function "*" (Left : Complex_Vector;

Right : Real_Vector) **return** Complex_Vector;

Each operation returns the vector product of Left and Right. The exception
Constraint_Error is raised if Left'Length and Right'Length are not equal to 3.
The index range of the result is Left'Range.

This function returns a "unit vector" with Order components and a lower
bound of First. All components are set to (0.0,0.0) except for the Index
component which is set to (1.0,0.0). The exception Constraint_Error is
raised if Index < First, Index > First + Order - 1, or
if First + Order - 1 > Integer'Last.

Each function returns a matrix of the specified cartesian components of X.
The index ranges of the result are those of X.

Each procedure replaces the specified (cartesian) component of each of the
components of X by the value of the matching component of Re or Im; the other
(cartesian) component of each of the components is unchanged. The exception
Constraint_Error is raised if X'Length(1) is not equal to Re'Length(1) or
X'Length(2) is not equal to Re'Length(2) and if X'Length(1) is not equal to
Im'Length(1) or X'Length(2) is not equal to Im'Length(2).

Each function constructs a matrix of Complex results (in cartesian
representation) formed from given matrices of cartesian components; when
only the real components are given, imaginary components of zero are
assumed. The index ranges of the result are those of Re. The exception
Constraint_Error is raised if Re'Length(1) is not equal to Im'Length(1) or
Re'Length(2) is not equal to Im'Length(2).

Each function calculates and returns a matrix of the specified polar
components of X or Right using the corresponding function in
Numerics.Generic_Complex_Types. The index ranges of the result are those of X
or Right.

Each function constructs a matrix of Complex results (in cartesian
representation) formed from given matrices of polar components using
the corresponding function in Numerics.Generic_Complex_Types on matching
components of Modulus and Argument. The index ranges of the result are those
of Modulus. The exception Constraint_Error is raised if Modulus'Length(1) is
not equal to Argument'Length(1) or Modulus'Length(2) is not equal to
Argument'Length(2).

function "+" (Right : Complex_Matrix) return Complex_Matrix;
function "-" (Right : Complex_Matrix) return Complex_Matrix;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of Right. The index
ranges of the result are those of Right.

This function returns the result of applying the appropriate function Conjugate
in Numerics.Generic_Complex_Types to each component of X. The index ranges of
the result are those of X.

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 : Complex_Matrix) return Complex_Matrix;
function "-" (Left, Right : Complex_Matrix) return Complex_Matrix;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of Left and the matching
component of Right. The index ranges of the result are those of
Left. The exception 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 : Complex_Matrix) return Complex_Matrix;

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. The exception
Constraint_Error is raised if Left'Length(2) is not equal to Right'Length(1).
This operation involves inner products.

function "*" (Left : Complex_Vector;

Right : Complex_Matrix) **return** Complex_Vector;

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). The exception Constraint_Error is raised if
Left'Length is not equal to Right'Length(1). This operation involves inner
products.

function "*" (Left : Complex_Matrix;

Right : Complex_Vector) **return** Complex_Vector;

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). The exception Constraint_Error is raised if
Left'Length(2) is not equal to Right'Length. This operation involves inner
products.

function "+" (Left : Real_Matrix;

Right : Complex_Matrix) **return** Complex_Matrix;

function "+" (Left : Complex_Matrix;

Right : Real_Matrix) **return** Complex_Matrix;

function "-" (Left : Real_Matrix;

Right : Complex_Matrix) **return** Complex_Matrix;

function "-" (Left : Complex_Matrix;

Right : Real_Matrix) **return** Complex_Matrix;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of Left and the matching
component of Right. The index ranges of the result are those of Left. The
exception 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 : Real_Matrix;

Right : Complex_Matrix) **return** Complex_Matrix;

function "*" (Left : Complex_Matrix;

Right : Real_Matrix) **return** Complex_Matrix;

Each 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. The exception
Constraint_Error is raised if Left'Length(2) is not equal to Right'Length(1).
These operations involve inner products.

function "*" (Left : Real_Vector;

Right : Complex_Matrix) **return** Complex_Vector;

function "*" (Left : Complex_Vector;

Right : Real_Matrix) **return** Complex_Vector;

Each 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). The exception Constraint_Error is raised if
Left'Length is not equal to Right'Length(1). These operations involve inner
products.

function "*" (Left : Real_Matrix;

Right : Complex_Vector) **return** Complex_Vector;

function "*" (Left : Complex_Matrix;

Right : Real_Vector) **return** Complex_Vector;

Each 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). The exception Constraint_Error is raised if
Left'Length(2) is not equal to Right'Length. These operations involve inner
products.

function "*" (Left : Complex; Right : Complex_Matrix) return Complex_Matrix;

This operation returns the result of multiplying each component of Right by
the complex number Left using the appropriate operation "*" in
Numerics.Generic_Complex_Types. The index ranges of the result are those of
Right.

function "*" (Left : Complex_Matrix; Right : Complex) return Complex_Matrix;
function "/" (Left : Complex_Matrix; Right : Complex) return Complex_Matrix;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of the matrix Left and the
complex number Right. The index ranges of the result are those of Left.

function "*" (Left : Real'Base; Right : Complex_Matrix) return Complex_Matrix;

This operation returns the result of multiplying each component of Right by
the real number Left using the appropriate operation "*" in
Numerics.Generic_Complex_Types. The index ranges of the result are those of
Right.

function "*" (Left : Complex_Matrix; Right : Real'Base) return Complex_Matrix;
function "/" (Left : Complex_Matrix; Right : Real'Base) return Complex_Matrix;

Each operation returns the result of applying the corresponding operation in
Numerics.Generic_Complex_Types to each component of the matrix Left and the
scalar Right. The index ranges of the result are those of Left.

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 X'Range. 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.

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 those of X. 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.

This function returns a matrix B such that A * B is (nearly) the unit matrix.
The index ranges of the result are those of A. 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.

This function returns the determinant of the matrix A. Constraint_Error is
raised if A'Length(1) is not equal to A'Length(2).

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,0.0) except for the main diagonal,
whose components are set to (1.0,0.0). The exception Constraint_Error is raised
if First_1 + Order - 1 > Integer'Last or First_2 + Order - 1 > Integer'Last.

Implementation Requirements

Accuracy requirements for the functions Solve, Inverse and Determinant are
implementation defined.

For operations not involving an inner product, the accuracy requirements are
those of the corresponding operations of the type Real'Base and Complex in both
the strict mode and the relaxed mode (see G.2).

?? This doesn't really cover vector product but I might delete it anyway.

For operations involving an inner product, no requirements are specified in
the relaxed mode. In the strict mode the acuracy should be at least that of
the canonical implementation of multiplication and addition using the
corresponding operations of type Real'Base and Complex and performing the
cumulative addition using ascending indices. Implementations shall document
any techniques used to reduce cancellation errors such as extended precision
arithmetic.

Implementation Permissions

The nongeneric equivalent packages may, but need not, be actual
instantiations of the generic package for the appropriate predefined type.

Although many operations are defined in terms of operations from
Numerics.Generic_Complex_Types, they need not be implemented by calling those
operations provided that the effect is the same.

Implementation Advice

Implementations should implement the Solve and Inverse functions using
established techniques such as Gauss-Seidel with Interchanges.

It is not the intention that any special provision should be made to
determine whether a matrix is ill-conditioned or not. The naturally occuring
division by zero or overflow which will result from executing these
functions with an ill-conditioned matrix and thus raise Constraint_Error is
sufficient.

Implementations should not perform operations on mixed complex and real operands
by first converting the real operand to complex. See G.1.1(56,57).

G.3.3 Eigensystems of Matrices

Eigenvalues and eigenvectors of a real symmetric matrix are computed using
subprograms in the child package Ada.Numerics.Generic_Real_Arrays.Eigen
whose declaration is:

The library package Numerics.Real_Arrays.Eigen is declared pure and defines
the same subprograms as Numerics.Generic_Real_Arrays.Eigen, 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.Eigen, Numerics.Long_Real_Arrays.Eigen, etc.

The effect of the subprograms is as described below.

This function returns the eigenvalues of the symmetric matrix A as a vector
sorted into order with the largest first. The exception 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).

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 Values. 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. The exception Constraint_Error is raised if A'Length(1) is not
equal to A'Length(2). The index ranges of the parameter Vectors are those of A.

For both subprograms, if the matrix A is not symmetric then Constraint_Error is
raised.

?? Or we could define the matrix just in terms of the upper diagonal or we
could make it symmetric by taking the arithmetic mean of the matrix and its
transpose.

Eigenvalues and eigenvectors of a general real or complex matrix are computed
using subprograms in the child package
Ada.Numerics.Generic_Complex_Arrays.Eigen whose declaration is:

-- Eigensystem of a general real matrix

-- Eigensystem of a Hermitian matrix

-- Eigensystem of a general complex matrix

The library package Numerics.Complex_Arrays.Eigen is declared pure and defines
the same subprograms as Numerics.Generic_Complex_Arrays.Eigen, except that
the Real_Vector and Real_Matrix types exported by Numerics.Real_Arrays are
systematically substituted for Real_Vector and Real_Matrix throughout.
Nongeneric equivalents for each of the other predefined floating point types
are defined similarly, with the names
Numerics.Short_Complex_Arrays.Eigen, Numerics.Long_Complex_Arrays.Eigen, etc.

?? I didn't mention Real'Base, Complex etc because they are not used in the
spec. Should I??

The effect of the subprograms is as described below.

This function returns the eigenvalues of the matrix A as a vector sorted by
order of modulus with the largest first. The exception 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).

This procedure computes both the eigenvalues and eigenvectors of the
matrix A. The out parameter Values is the same as that obtained by calling the
function Values. 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 exception Constraint_Error is
raised if A'Length(1) is not equal to A'Length(2). The index ranges of the
parameter Vectors are those of A.

This function returns the eigenvalues of the hermitian matrix A as a vector
sorted into order with the largest first. The exception 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). The exception Constraint_Error is raised if the matrix
is not hermitian.

This procedure computes both the eigenvalues and eigenvectors of the hermitian
matrix A. The out parameter Values is the same as that obtained by calling the
function Values. 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 mutually orthonormal,
including when there are repeated eigenvalues. The exception Constraint_Error
is raised if A'Length(1) is not equal to A'Length(2). The index ranges of the
parameter Vectors are those of A. The exception Constraint_Error is raised if
the matrix is not hermitian.

?? Or we could define the matrix just in terms of the upper diagonal or we
could make it hermitian by taking the arithmetic mean of the matrix and its
conjugate transpose.

This function returns the eigenvalues of the matrix A as a vector sorted by
order of modulus with the largest first. The exception 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).

This procedure computes both the eigenvalues and eigenvectors of the
matrix A. The out parameter Values is the same as that obtained by calling the
function Values. 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 exception Constraint_Error is
raised if A'Length(1) is not equal to A'Length(2). The index ranges of the
parameter Vectors are those of A.

Implementation Requirements

The accuracy of these subprograms is implementation defined.

Implementation Permissions

The nongeneric equivalent packages may, but need not, be actual
instantiations of the generic package for the appropriate predefined type.

!example

To be done by someone else.

!discussion

The packages closely follow those in 13813. However, the discussion has been
considerably simplified by assuming the properties of the corresponding scalar
operations such as those in Numerics.Complex_Types. This removes a lot of
explicit mention of raising exceptions and accuracy. Also remarks that these
are standard mathematical operations have been deleted when the meaning is
given by other words in the description.

The component by component operations of * / and ** on vectors have been
deleted on the grounds that they are not useful. (They might be useful for
manipulating arrays in general but we are concerned with arrays used as vectors
for linear algebra.)

Operations for vector products have been added, However, this might be
considered curious because they only apply in three-dimensional space.

?? And I put the inner outer and vector products together which might not have
been a good idea.

It is hoped that there is not too much confusion between component when
applied to the parts of a complex number and component when applied to a part
of an array. 13813 uses element for the latter but the proper Ada term for an
element of an array or field of a record is of course component.

Observe that the index range of the various arrays is Integer (rather than
Natural or Positive). This permits negative indices and in particular arrays
with symmetric index ranges about zero.

The function Identity_Matrix of 13813 has been changed to Unit_Matrix on the
grounds that the prime concern is with linear algebra and not group theory.

Functions have been added to Numerics.Generic_Real_Arrays for matrix inversion
and related operations. On the grounds that it would be helpful to have simple
algorithms available to users rather than none at all, no accuracy
requirements are specified. Instead the accuracy is stated to be
implementation-defined which means that the implementation must document it.

The names chosen are Solve and Inverse. The former is overloaded, one version
solves a single set of linear equations; the second solves multiple sets. Note
that Inverse is not strictly necesssary becuase the effect can be obtained by
using Solve and a Unit_Matrix thus

I := Unit_Matrix(A'Length); B := Solve(A, I); -- same as B := Inverse (A);

A function for computing the determinant has been added since it is helpful in
deciding whether an array is ill-conditioned and therefore the results of
inversion might be suspect.

Similar functions have also been added for complex arrays.

In addition, subprograms have been added for the computation of eigenvalues
and vectors of real and complex arrays. Because these are somewhat more
specialized than the elementary operations they have been placed in generic
child packages such as Numerics.Generic_Real_Arrays.Eigen. The individual
sunprograms are Values and Values_And_Vectors which means that they can be
invoked by statements such as

V := Eigen.Values(A);

There is no separate function for eigenvectors since it is unlikely that these
would be required without the eigenvalues.

The most common application is with real symmetric matrices and these have real
eigenvalues and vectors. Hence the subprograms in the real package only cater
for symmetric matrices. Applications include: moments of inertia, elasticity,
confidence regions and so on.

A general real matrix can have complex eigenvalues and therefore provision for
these and complex matrices is in the complex package. Note that the analogue of
symmetric matrices in the complex domain is hermitian matrices (the transpose
of a hermitian matrix equals its complex conjugate); these have real
eigenvalues and distinct provision is made for these. The various overloadings
are conveniently distinguished by the parameter profiles. (but could not be if
there were separate fuinctions for the eigenvectors).

?? I am somewhat rusty on eigensystems of non-symmetric matrices. I have said
nothing about the normalization of the vectors. Note that I put all the eigen
stuff in 3.3 - that's partly so it can be removed easily!!

The original text of 13813 contains the remark "no complex conjugation is
performed" in the description of multiplication of complex vectors to produce
a complex matrix. I removed it because it seemed unhelpful to the user. The
origin of this according to Donald Sando is that the NUMWG considered including
various combined operations involving both multiplication and conjugation and
indeed multiplication and transposition. Multiplication and conjugation are
natural when dealing with for example Hermitian matrices. however, to provided
all the appropriate combinations would have greatly increased the size of the
package so it says in the rationale section F.3 of the 1998 version of 13813.

Terry Froggatt adds further thoughts on this matter (private communication,
17 July 2001. Thus

2. THOUGHTS ON TRANSPOSITION.

This comment arises from the last paragraph of F.3 of the Rationale in
Annexe F ....

In the Harrier flight-program we have several spatial co-ordinate systems,
for example:

Earth Axes: North - East - Up.
Heading - Cross-heading - Up, (yawed wrt Earth).
Aircraft Axes (rolled/pitched/yawed wrt Earth).
HUD Axes (6 degrees down from Aircraft Axes).

Conversion of vectors between pairs of these systems uses matrix
multiplication one way, and multiplication by the tranpose of the same
matrix the other way.

Whilst X*CONJUGATE(Y) may be acceptable, since X and Y are the same size,
the costs for X*TRANSPOSE(Y) are different: Y is generally much bigger than
X, and it seems rather wasteful (mainly in time) to put the transpose of a
matrix onto the stack then throw it away.

An alternative is to calculate T := TRANSPOSE(Y) as soon as Y is calculated,
so that it is available to convert several vectors between a pair of axis
systems. This is wasteful (mainly in space) because several transposes have
to be kept around simultaneously.

What we actually use is one routine which implements X*TRANSPOSE(Y) directly:
it performs the matrix multiplication but traverses Y in transpose-order. This
takes no more time than the normal multiply, and no space for transposes (just
code space).

This is exactly what 13813's rationale rejected.

Note that it should be possible to express the Harrier routine, and calls to
it, quite neatly in Ada. Following the style of "Imaginary", define:

function "*" (Left: Real_Vector; Right: Transposed_Real_Matrix)

return (Real_Vector);

This "*" knows it has to traverse Right in transposed order, and can be
called (using operator notation, and no qualification) by

"Z := X*Transpose_In_Situ(Y);".

Well now. It seems to me that playing with Hermitian matrices is pretty
rare and so we can forget about possible combinations of conjugate and
multiply. However, maybe there is a case for adding some combined
multiplication and transposition for the real matrices as explained above.

!ACATS Test

ACATS test(s) need to be created.

!appendix

**************************************************************** Recommendation on ISO/IEC 13813 from the UK The standard ISO/IEC 13813 entitled generic packages of real and complex type declarations and basic operations for Ada (including vector and matrix types) will soon be up for review. This note reviews the background to the development of the standard and makes a recommendation that the standard be revised. Background The Numerics Working Group of WG9 met many times during the period when Ada 95 was being designed and produced a number of standards. They were faced with the problem of whether to produce standards based on Ada 83 (87 in ISO terms) or whether to base them on Ada 95 or subsume them into Ada 95. One dilemma was of course that although Ada 95 was on the way nevertheless Ada 83 was expected to continue in use for many years. The standards are 11430: Generic package of elementary functions for Ada. 11729: Generic package of primitive functions for Ada. 13813: Generic packages 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. 11430 and 11729 are mentioned for completeness. They were published in 1994. They were based entirely on Ada 83 and their facilities are provided in the Ada 95 core language. The elementary functions, 11430, became the package Ada.Numerics.Generic_Elementary_Functions and the primitive functions, 11729, became the various attributes such as 'Floor and 'Ceiling, and 'Exponent and 'Fraction. These two standards were withdrawn recently and need no further mention. The other two standards, 13813 and 13814, were published in 1998 and will soon be up for review at the end of their five year period. Three possible fates can befall a standard when it is reviewed. It can be withdrawn, revised or confirmed. In the case of 13814, the functionality is all incorporated into the Numerics Annex of Ada 95 as the package Ada.Numerics.Generic_Complex_Elementary_Functions. There are a few changes in presentation because the Ada 95 package uses the generic package parameter feature which of course did not exist in Ada 83. Nevertheless there seems little point in continuing with 13814 and so at the Leuven meeting of WG9 it was agreed to recommend that it be withdrawn. However, the situation regarding 13813 is not so clear. Some of its functionality is included in Ada 95 but quite a lot is not. The topics covered are (1) a complex types package including various complex arithmetic operations, (2) a real arrays package covering both vectors and matrices, (3) a complex arrays package covering both vectors and matrices, (4) a complex input-output package. The complex types package (1) became the package Ada.Numerics.Generic_Complex_Types and the input-output package (4) became Ada.Text_IO.Complex_IO. However, the array packages, both real and complex, were not incorporated into the Ada 95 standard. At the Leuven meeting, it was agreed that 13813 should not be withdrawn without further study. The UK was asked to study whether small or large changes are required in 13813 and to report back. The Ada Rapporteur Group would then decide whether the functionality should be included in a future revision or amendment to Ada 95. This is the report from the UK. Recommendation It is recommended that 13813 be revised so that it only contains the functionality not included in Ada 95. The revised standard should contain two generic packages namely Ada.Numerics.Generic_Real_Arrays and Ada.Numerics.Generic_Complex_Arrays. There should also be standard non-generic packages corresponding to the predefined types such as Float in an analogous manner to the standard packages such as Ada.Numerics.Complex_Types and Ada.Numerics.Long_Complex_Types for Float and Long_Float respectively. The text of the Ada specifications of the two generic packages should be essentially as given in the nonnormative Annex G of the existing standard 13813. (This Annex illustrates how the existing standard packages might be rewritten using Ada 95. There is an error regarding the formal package parameters which has been corrected in the revised text.) There is an important issue regarding what should happen if there is a mismatch in the array lengths of the parameters in a number of the subprograms provided by the packages. For example if function "+" (Left, Right: Real_Vector) return Real_Vector be called with parameters such that Left'Length /= Right'Length. The existing standard raises the exception Array_Index_Error which is declared alone in a package Array_Exceptions. The nonnormative Annex G shows this exception incorporated into the package Ada.Numerics thereby producing an incompatibility with the existing definition of Ada.Numerics. We considered four possibilities regarding this exception 1) Add Array_Index_Error to Ada.Numerics as in Annex G. 2) Place Array_Index_Error in a new child package such as Ada.Numerics.Arrays. 3) Eliminate Array_Index_Error and raise Constraint_Error instead. 4) State that the behaviour with mismatched arrays is implementation-defined. We concluded that Option (1) is undesirable because of incompatibility. Option (2) is feasible but one ought then to place the generic packages themselves into this package so that they become Ada.Numerics.Arrays.Generic_Real_Arrays and Ada.Numerics.Arrays.Generic_Complex_Arrays. This nesting is considered cumbersome. Option(3) gives the same behaviour as similar mismatching on predefined operations and although losing some specificacity has practical simplicity. Option (4) is disliked since gratuitous implementation-defined behaviour should be avoided. We therefore recommend Option (3) that Constraint_Error be raised on mismatching of parameters. If the Ada95 standard itself be revised at some later date then consideration should be given to incorporating the functionality of the revised 13813 into the Numerics Annex. Proposed text The proposed normative text of the revised standard is distributed separately as N404, Working Draft, Revision of ISO/IEC 13813. Note that the non-normative rationale section remains to be completed and that consequently the bibliography might need alterations to match. Acknowledgment We acknowledge the valuable assistance of Donald Sando, the editor of the original standard, in the preparation of this recommendation. **************************************************************** From: John Barnes Sent: Tuesday, January 21, 2003 10:33 AM Gosh - I have made a start on AI-296. [Editor's note: This is version /01 of the AI.] I have pulled all the Ada text in and written a first cut at the description for the real vectors and matrices. I thought maybe it would be a good idea to get the style settled before spending time on the complex ones. Here are some thoughts. There are questions of how to arrange the annex in !proposal. What do we use for the plural of index when talking about arrays? I suspect that it ought to be indexes and not indices as Don had. I have slimmed down the accuracy and error stuff that Don Sando had by simply referring to the underlying real operations. There are as yet no useful remarks on the accuracy of inner product. I started by explaining the bounds of the result in terms of the bounds of the parameters whereas Don had done it in terms of ranges. Later I discovered that ranges are much easier for the more elaborate matrix cases so perhaps I should have used ranges throughout. Thus saying for example "the index range of the result is Left'Range" rather than "the bounds of the result are those of the left operand". I am not sure whether I need to say anything about null arrays. I note that 4.5.1(8) doesn't seem to. Incidentally, I am still overwhelmed with other work. Although the Spark book has gone away to be proof read, it re-emerges tomorrow and will keep me busy for several days in doing final corrections. Also I have to prepare a one-day course on Spark at the University of York the week after Padua so I have to send the notes beforehand. I do this course each year but at the last minute they want a lot of changes. But we do have a baseline to discuss on AI-296 for Padua even if I don't get a chance to spend much more time on it before then. So I don't feel too guilty. There is a lot of interest in this stuff (especially problems of accuracy) in the UK and we have a BSI meeting in mid February to discuss this. **************************************************************** From: Pascal Leroy Sent: Wednesday, January 22, 2003 8:24 AM > There are as yet no useful remarks on the > accuracy of inner product. I think that the accuracy of the inner product should be defined in a way similar to that of the complex multiplication. If you look at the inner product of vectors (a1, a2) and (b1, b2), the result is quite similar to the real part of the (complex) multiplication of a1 + ib1 and a1 - ib2. The box error in G.2.5 ensures that if the result of a multiplication lands very close to one of the axis, you don't have to provide an unreasonable accuracy on the "small" component, because presumably cancellations happened. Moreover, the "large" component gives you some information on the magnitude of the cancellation, so it can be used to define the acceptable error. That's essentially what the box error does (at least that's my understanding). For the inner product of v1 and v2, I believe that the error should similarly be defined to be of the form d * length(v1) * length (v2). If cancellations happen, i.e. if the vectors are nearly orthogonal, the error can be quite large compared to the final result. In the absence of cancellations, i.e. if the vectors are not orthogonal, then with a proper choice for d we can require good accuracy. The alert reader will notice that this amounts to defining the relative error on the result as d/cos (A) where A is the angle of the two vectors. This is just a back-of-an-envelope proposal. No rigorous analysis was done. **************************************************************** From: John Barnes Sent: Tuesday, June 10, 2003 4:42 AM Vectors and matrices At the last ARG meeting I was asked to add material for matrix inversion, solution of linear equations, determinants, and eigenvalues and vectors. We very briefly considered a couple of other topics but dismissed them as tricky. I have updated AI-296 as instructed and sent it to Randy so I assume it is now on the database. [Editor's note: This is version /02 of the AI.] I put the linear equations, invert and determinant in the main packages (one for reals and one for complex). I put the eigenvalues and vectors in child packages because they seemed a bit different - also they were more complex than I had remembered (I hadn't looked at this stuff for nearly half a century) and maybe they are more elaborate than we need. There are a number of queries embedded in the text concerning accuracy etc. However, there are a number of other matters which ought to be discussed. I visited NAG (National Algorithms Group). They did a lot of Ada 83 numerics stuff once although they no longer sell it. They did implement 13813. They put the invert stuff in separate packages but it is much more elaborate than my humble proposal. For example they have options on whether to invert just once or to iterate on the residuals until they are acceptable. But then of course one has to give the accuracy required (although they do have defaults). Moreover the LU decomposition is made visible. Now as a user, I don't really want to know about how it is done (in a sense it breaks the abstraction), I just want the answer. They have provided what the numerical analyst specialist would want and I don't really think that is appropriate for a simple standard but maybe I am wrong. Thus their subprograms have some 10 parameters whereas in my proposal they are the minimal one or two. For instance their determinant has eight parameters whereas mine has one. In the case of eigenvalues and vectors, they had options to select just the largest or smallest or all those in a certain range. No doubt because of timing considerations for large matrices. We also considered other topics which might provide candidates for a standard. The trouble with many topics is that techniques are still evolving and it is unclear that a unique worthy proposal would be easy. Curve and surface fitting. There is lots of choice here. Cubic splines are fashionable. Chebychev polynolmials are a long established technique. But both require significant user knowledge and it is easy to misuse them. Maximization/minimization. Lots of options here and real problems have constraints. The techniques are heuristic and still evolving so this looks really hard and can be eliminated. Partial differential equations seem too much for the casual user. In any event there is a huge variety of parameters and constraints. Ordinary differential equations and quadrature (integration) are perhaps more promising. Again there are quite a lot of options. But this might be worth exploring. Roots of general non-linear equations. An evolving subject still. Very heuristic and prone to chaotic difficulties. Roots of polynomial equations. This seemed more likely . It appears that the subject is static and a well-established technique devised by Laguerre is typically used. (But it is outside my domain of knowledge.) I have also contemplated my 20 year old HP pocket calculator. I had a vague feeling that Ada should be able to do what it can do. It can do the following: The usual trigonometric and hyperbolic functions. Real and complex matrices, linear equations, inversion and determinants (reveals the LU decomposition). It does not do eigenvalues. Random number generation. A whole heap of statistical stuff. Integration (quadrature) of a function of one variable. Solution of a non-linear equation in one variable. So where do we go from here? It is a slippery slope and I fear that I am not an expert. ****************************************************************

Questions? Ask the ACAA Technical Agent