!standard G.3 (01) 04-05-24 AI95-00296/09
!standard G.3.1 (01)
!standard G.3.2 (01)
!class amendment 02-06-07
!status Amendment 200Y 04-01-13
!status ARG Approved 12-0-1 03-12-11
!status work item 03-10-29
!status ARG Approved 10-0-0 03-10-03
!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
subprograms for the solution of linear equations, matrix inversion,
determinants, and the determination of the eigenvalues and eigenvectors of
real symmetric matrices and Hermitian matrices.
!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:
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;
-- 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;
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. Constraint_Error
is raised if Left'Length is not equal to Right'Length.
function "*" (Left, Right : Real_Vector) return Real'Base;
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 "*" (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 Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) return Real_Vector;
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;
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;
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
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;
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;
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
matrix result are Left'Range and Right'Range respectively.
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). 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). 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.
function Solve (A : Real_Matrix; X: Real_Vector) return Real_Vector;
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.
function Solve (A, X : Real_Matrix) return Real_Matrix;
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.
function Inverse (A : Real_Matrix) return Real_Matrix;
This function returns a matrix B such that A * B is (nearly) equal to 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.
function Determinant (A : Real_Matrix) return Real'Base;
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;
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);
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). The index ranges of the parameter Vectors are those of A.
Argument_Error is raised if the matrix A is not symmetric.
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1) return Real_Matrix;
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
Accuracy requirements for the subprograms Solve, Inverse, Determinant,
Eigenvalues and Eigensystem 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).
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'Machine_Mantissa)
Documentation Requirements
Implementations shall document any techniques used to reduce cancellation
errors such as extended precision arithmetic.
AARM 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.
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 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.
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.
The test that a matrix is symmetric may be performed by using the equality
operator to compare the relevant components.
G.3.2 Complex Vectors and Matrices
The generic library package Numerics.Generic_Complex_Arrays has the following
declaration:
with Ada.Numerics.Generic_Real_Arrays, Ada.Numerics.Generic_Complex_Types;
generic
with package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (<>);
use Real_Arrays;
with package Complex_Types is new Ada.Numerics.Generic_Complex_Types (Real);
use Complex_Types;
package Ada.Numerics.Generic_Complex_Arrays is
pragma Pure(Generic_Complex_Arrays);
-- Types
type Complex_Vector is array (Integer range <>) of Complex;
type Complex_Matrix is array (Integer range <>,
Integer range <>) of Complex;
-- Subprograms for Complex_Vector types
-- Complex_Vector selection, conversion and composition operations
function Re (X : Complex_Vector) return Real_Vector;
function Im (X : Complex_Vector) return Real_Vector;
procedure Set_Re (X : in out Complex_Vector;
Re : in Real_Vector);
procedure Set_Im (X : in out Complex_Vector;
Im : in Real_Vector);
function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector;
function Compose_From_Cartesian (Re, Im : Real_Vector) return Complex_Vector;
function Modulus (X : Complex_Vector) return Real_Vector;
function "abs" (Right : Complex_Vector) return Real_Vector
renames Modulus;
function Argument (X : Complex_Vector) return Real_Vector;
function Argument (X : Complex_Vector;
Cycle : Real'Base) return Real_Vector;
function Compose_From_Polar (Modulus, Argument : Real_Vector)
return Complex_Vector;
function Compose_From_Polar (Modulus, Argument : Real_Vector;
Cycle : Real'Base)
return Complex_Vector;
-- Complex_Vector arithmetic operations
function "+" (Right : Complex_Vector) return Complex_Vector;
function "-" (Right : Complex_Vector) return Complex_Vector;
function Conjugate (X : Complex_Vector) return Complex_Vector;
function "+" (Left, Right : Complex_Vector) return Complex_Vector;
function "-" (Left, Right : Complex_Vector) return Complex_Vector;
function "*" (Left, Right : Complex_Vector) return Complex;
-- Mixed Real_Vector and Complex_Vector arithmetic operations
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;
function "*" (Left : Real_Vector; Right : Complex_Vector) return Complex;
function "*" (Left : Complex_Vector; Right : Real_Vector) return Complex;
-- Complex_Vector scaling operations
function "*" (Left : Complex;
Right : Complex_Vector) return Complex_Vector;
function "*" (Left : Complex_Vector;
Right : Complex) return Complex_Vector;
function "/" (Left : Complex_Vector;
Right : Complex) return Complex_Vector;
function "*" (Left : Real'Base;
Right : Complex_Vector) return Complex_Vector;
function "*" (Left : Complex_Vector;
Right : Real'Base) return Complex_Vector;
function "/" (Left : Complex_Vector;
Right : Real'Base) return Complex_Vector;
-- Other Complex_Vector operations
function Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) return Complex_Vector;
-- Subprograms for Complex_Matrix types
-- Complex_Matrix selection, conversion and composition operations
function Re (X : Complex_Matrix) return Real_Matrix;
function Im (X : Complex_Matrix) return Real_Matrix;
procedure Set_Re (X : in out Complex_Matrix;
Re : in Real_Matrix);
procedure Set_Im (X : in out Complex_Matrix;
Im : in Real_Matrix);
function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix;
function Compose_From_Cartesian (Re, Im : Real_Matrix) return Complex_Matrix;
function Modulus (X : Complex_Matrix) return Real_Matrix;
function "abs" (Right : Complex_Matrix) return Real_Matrix
renames Modulus;
function Argument (X : Complex_Matrix) return Real_Matrix;
function Argument (X : Complex_Matrix;
Cycle : Real'Base) return Real_Matrix;
function Compose_From_Polar (Modulus, Argument : Real_Matrix)
return Complex_Matrix;
function Compose_From_Polar (Modulus, Argument : Real_Matrix;
Cycle : Real'Base)
return Complex_Matrix;
-- Complex_Matrix arithmetic operations
function "+" (Right : Complex_Matrix) return Complex_Matrix;
function "-" (Right : Complex_Matrix) return Complex_Matrix;
function Conjugate (X : Complex_Matrix) return Complex_Matrix;
function Transpose (X : Complex_Matrix) return Complex_Matrix;
function "+" (Left, Right : Complex_Matrix) return Complex_Matrix;
function "-" (Left, Right : Complex_Matrix) return Complex_Matrix;
function "*" (Left, Right : Complex_Matrix) return Complex_Matrix;
function "*" (Left, Right : Complex_Vector) return Complex_Matrix;
function "*" (Left : Complex_Vector;
Right : Complex_Matrix) return Complex_Vector;
function "*" (Left : Complex_Matrix;
Right : Complex_Vector) return Complex_Vector;
-- Mixed Real_Matrix and Complex_Matrix arithmetic operations
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;
function "*" (Left : Real_Matrix;
Right : Complex_Matrix) return Complex_Matrix;
function "*" (Left : Complex_Matrix;
Right : Real_Matrix) return Complex_Matrix;
function "*" (Left : Real_Vector;
Right : Complex_Vector) return Complex_Matrix;
function "*" (Left : Complex_Vector;
Right : Real_Vector) return Complex_Matrix;
function "*" (Left : Real_Vector;
Right : Complex_Matrix) return Complex_Vector;
function "*" (Left : Complex_Vector;
Right : Real_Matrix) return Complex_Vector;
function "*" (Left : Real_Matrix;
Right : Complex_Vector) return Complex_Vector;
function "*" (Left : Complex_Matrix;
Right : Real_Vector) return Complex_Vector;
-- Complex_Matrix scaling operations
function "*" (Left : Complex;
Right : Complex_Matrix) return Complex_Matrix;
function "*" (Left : Complex_Matrix;
Right : Complex) return Complex_Matrix;
function "/" (Left : Complex_Matrix;
Right : Complex) return Complex_Matrix;
function "*" (Left : Real'Base;
Right : Complex_Matrix) return Complex_Matrix;
function "*" (Left : Complex_Matrix;
Right : Real'Base) return Complex_Matrix;
function "/" (Left : Complex_Matrix;
Right : Real'Base) return Complex_Matrix;
-- Complex_Matrix inversion and related operations
function Solve (A : Complex_Matrix; X: Complex_Vector) return Complex_Vector;
function Solve (A, X : Complex_Matrix) return Complex_Matrix;
function Inverse (A : Complex_Matrix) return Complex_Matrix;
function Determinant (A : Complex_Matrix) return Complex;
-- Eigenvalues and vectors of a Hermitian matrix
function Eigenvalues (A : Complex_Matrix) return Real_Vector;
procedure Eigensystem (A : in Complex_Matrix;
Values : out Real_Vector;
Vectors : out Complex_Matrix);
-- Other Complex_Matrix operations
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1)
return Complex_Matrix;
end Ada.Numerics.Generic_Complex_Arrays;
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 type exported by Numerics.Complex_Types is systematically
substituted for Complex, 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.
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 are 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.
function Re (X : Complex_Vector) return Real_Vector;
function Im (X : Complex_Vector) return Real_Vector;
Each function returns a vector of the specified cartesian components of X.
The index range of the result is X'Range.
procedure Set_Re (X : in out Complex_Vector; Re : in Real_Vector);
procedure Set_Im (X : in out Complex_Vector; Im : in Real_Vector);
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.
Constraint_Error is raised if X'Length is not equal to
Re'Length or Im'Length.
function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector;
function Compose_From_Cartesian (Re, Im : Real_Vector) return Complex_Vector;
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. Constraint_Error is
raised if Re'Length is not equal to Im'Length.
function Modulus (X : Complex_Vector) return Real_Vector;
function "abs" (Right : Complex_Vector) return Real_Vector renames Modulus;
function Argument (X : Complex_Vector) return Real_Vector;
function Argument (X : Complex_Vector;
Cycle : Real'Base) return Real_Vector;
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.
function Compose_From_Polar (Modulus, Argument : Real_Vector)
return Complex_Vector;
function Compose_From_Polar (Modulus, Argument : Real_Vector; Cycle : Real'Base)
return Complex_Vector;
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. 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.
function Conjugate (X : Complex_Vector) return Complex_Vector;
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;
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.
Constraint_Error is raised if Left'Length is not equal to
Right'Length.
function "*" (Left, Right : Complex_Vector) return Complex;
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 "+" (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.
Constraint_Error is raised if Left'Length is not equal to
Right'Length.
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.
Constraint_Error is raised if Left'Length is not equal to Right'Length.
These operations involve an inner product.
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 Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) return Complex_Vector;
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). Constraint_Error is
raised if Index < First, Index > First + Order - 1, or
if First + Order - 1 > Integer'Last.
function Re (X : Complex_Matrix) return Real_Matrix;
function Im (X : Complex_Matrix) return Real_Matrix;
Each function returns a matrix of the specified cartesian components of X.
The index ranges of the result are those of X.
procedure Set_Re (X : in out Complex_Matrix; Re : in Real_Matrix);
procedure Set_Im (X : in out Complex_Matrix; Im : in Real_Matrix);
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.
Constraint_Error is raised if X'Length(1) is not equal to Re'Length(1) or
Im'Length(1) or if X'Length(2) is not equal to Re'Length(2) or Im'Length(2).
function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix;
function Compose_From_Cartesian (Re, Im : Real_Matrix) return Complex_Matrix;
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.
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).
function Modulus (X : Complex_Matrix) return Real_Matrix;
function "abs" (Right : Complex_Matrix) return Real_Matrix renames Modulus;
function Argument (X : Complex_Matrix) return Real_Matrix;
function Argument (X : Complex_Matrix;
Cycle : Real'Base) return Real_Matrix;
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.
function Compose_From_Polar (Modulus, Argument : Real_Matrix)
return Complex_Matrix;
function Compose_From_Polar (Modulus, Argument : Real_Matrix;
Cycle : Real'Base)
return Complex_Matrix;
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. 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.
function Conjugate (X : Complex_Matrix) return Complex_Matrix;
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.
function Transpose (X : Complex_Matrix) return Complex_Matrix;
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. 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.
Constraint_Error is raised if Left'Length(2) is not equal to Right'Length(1).
This operation involves inner products.
function "*" (Left, Right : Complex_Vector) return Complex_Matrix;
This 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 : 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). 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). 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.
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_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 : 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). 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). 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
real number Right. The index ranges of the result are those of Left.
function Solve (A : Complex_Matrix; X: Complex_Vector) return Complex_Vector;
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.
function Solve (A, X : Complex_Matrix) return Complex_Matrix;
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.
function Inverse (A : Complex_Matrix) return Complex_Matrix;
This function returns a matrix B such that A * B is (nearly) equal to 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.
function Determinant (A : Complex_Matrix) return Complex;
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 : Complex_Matrix) return Real_Vector;
This function returns the eigenvalues of the Hermitian 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 Hermitian.
procedure Eigensystem (A : in Complex_Matrix;
Values : out Real_Vector;
Vectors : out Complex_Matrix);
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 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 mutually orthonormal,
including when there are repeated eigenvalues. 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. Argument_Error is raised if the
matrix A is not Hermitian.
function Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1)
return Complex_Matrix;
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). Constraint_Error is raised
if First_1 + Order - 1 > Integer'Last or First_2 + Order - 1 > Integer'Last.
Implementation Requirements
Accuracy requirements for the subprograms Solve, Inverse, Determinant,
Eigenvalues and Eigensystem 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).
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'Machine_Mantissa) for mixed
complex and real operands
g = sqrt(2.0) * X'Length * Real'Machine_Radix**(1-Real'Machine_Mantissa) for
two complex operands
Documentation Requirements
Implementations shall document any techniques used to reduce cancellation
errors such as extended precision arithmetic.
AARM 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 Complex 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),
Sections 3.1 and 3.6.
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. Implementations are recommended to refine the result by
performing an iteration on the residuals; if this is done then it should be
documented.
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.
The test that a matrix is Hermitian may use the equality operator to compare
the real components and negation followed by equality to compare the imaginary
components (see G.2.1).
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).
!example
The packages are self-explanatory and so no example is provided.
!discussion
Section G.1.1 of the Rationale for Ada 95 says
"A decision was made to abbreviate the Ada 95 packages by omitting the vector
and matrix types and operations. One reason was that such types and operations
were largely self-evident, so that little real help would be provided by
defining them in the language. Another reason was that a future version of Ada
might add enhancements for array manipulation and so it would be inappropriate
to lock in such operations prematurely."
It is now clear that such enhancements will not be added so the second
justification for omitting the facilities of 13813 disappears. In order to
overcome the objection that everything is self-evident we have taken the
approach that we should add some basic facilities that are commonly required,
not completely trivial to implement but on the other hand are mathematically
well understood.
The overall goal is thus twofold
* to provide commonly required facilities for the user who is not a numerical
professional,
* to provide a baseline of types and operations that forms a firm foundation
for binding to more general facilities such as the well-known BLAS (Basic
Linear Algebra Subprograms, see www.netlib.org/blas).
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 were considered but not added. This is because,
as usually defined, they only apply in three-dimensional space.
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.
The accuracy of most simple operations follows from the underlying operations
on scalar components. In the case of inner product there is the potential for
special operations to improve the speed and/or accuracy. We have specified
reasonable requirements in the strict mode, which are met by the canonical
implementation using a loop statement. This is because on some
hardware, built-in instructions which are fast actually lose accuracy. Note
that the Fortran language standard recognizes the existence of inner product
as a special case but imposes no accuracy requirements at all.
These requirements for inner product are based on the analysis in Chapter 3
of Accuracy and Stability of Numerical Algorithms by N J Higham. A factor of
nearly 2 has been added as a precaution. The mixed real and complex case uses
the same formula as the pure real case but a further factor of sqrt(2) is
required in the pure complex case essentially because of the introduction of
cross terms.
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 necessary because 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 common technique when solving sets of linear equations is to refine the result
by an iteration on the residuals. Thus to solve the set of equations Ax = y we
first perform
X := Solve(A, Y);
we then compute the error D by multiplying back thus
D := Y - A*X;
and then compute the refinement to X
DX := Solve(A, D);
and then correct X thus
X := X + DX;
Implementations are recommended to do this automatically; it requires little
extra computation if LU decomposition is used internally. If they do such
refinement then it should be documented.
A function for computing the determinant has been added since it is of some
help in deciding whether an array is ill-conditioned and therefore the results
of inversion might be suspect. Determinants are also useful in some statistical
calculations. The evaluation of determinants is very liable to overflow and many
such routines return a scaling power of 10 in order to keep the basic result
within range. For simplicity, it was decided not to do this since it is less
likely with matrices of low order; the user can of course scale the components
of the matrix if necessary.
Similar functions have also been added for complex arrays. However, it was not
deemed necessary to provide for mixed real and complex operands for Solve.
In addition, subprograms have been added for the computation of eigenvalues
and vectors of real symmetric matrices and Hermitian matrices. The subprograms
are Eigenvalues and Eigensystem. There is no separate subprogram 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. Applications include: moments of inertia, elasticity,
confidence regions and so on.
A slight problem arises in deciding who should check that a matrix is symmetric,
the user or the package. Computations of, for example, a covariance matrix will
result in a matrix that ought to be exactly symmetric but small errors might be
introduced which mean that it is not exactly symmetric. The onus could be
placed on the user to ensure that the matrix is exactly symmetric or
alternatively some tolerance could be passed to the eigen subprograms so that
they can perform a reasonable check. Passing a tolerance level adds an
irritating parameter, raises the issue of how the tolerance should be expressed
and gives the user the problem of what tolerance to request. Moreover, the
algorithms still have to decide which actual values should be used for those
pairs of components that are not exactly equal.
Accordingly, we have placed the onus on the user to ensure that the matrix is
exactly symmetric. This could be done for example by taking the mean of the
matrix and its transpose. The test in the subprograms can then test for exact
equality which is guaranteed to deliver the correct answer. If the test fails
then Argument_Error is raised. This exception, declared in Ada.Numerics, is
the normal exception for numeric operations when the argument is out of range.
Note that errors such as when bounds of arrays do not match raise
Constraint_Error by analogy with built-in array operations.
A third approach considered was for the user to supply a parameter indicating
which half of the matrix should be used to define it (the upper or lower
triangle). This avoids the need for any testing but it was considered bad
practice to be able to pass junk in the other half of the matrix.
The complex equivalent of real symmetric matrices are Hermitian matrices.
Hermitian matrices are such that their transpose (that is with rows and columns
interchanged) equals their complex conjugate (that is with the sign of imaginary
parts reversed). Hermitian matrices also have real eigenvalues.
Applications include Quantum Mechanics.
Again we have placed the onus on the user to ensure that the matrix is
Hermitian. The check in the package can then use strict equality for the real
parts and negation followed by equality for the imaginary parts.
We considered providing subprograms for the determination of eigenvalues and
eigenvectors of general real and complex matrices. Such matrices can have
complex eigenvalues and therefore provision for these would have to be in the
complex package. However, there are mathematical difficulties with these general
cases which are in strong contrast to the real symmetric and Hermitian matrices.
Thus, Numerical Recipes by Press, Flannery, Teukolsky and Vetterling says
regarding the real case:
"The algorithms for symmetric matrices ... are highly satisfactory in practice.
By contrast, it is impossible to design equally satisfactory algorithms for the
nonsymmetric case. There are two reasons for this. First, the eigenvalues of a
nonsymmetric matrix can be very sensitive to small changes in the matrix
elements. Second, the matrix itself can be defective so that there is no
complete set of eigenvectors. We emphasize that these difficulties are intrinsic
properties of certain nonsymmetric matrices, and no numerical procedure can cure
them."
Similar remarks apply to complex matrices where Hermitian matrices are well-
behaved but non-Hermitian matrices can be troublesome.
In view of these computational difficulties and the fact that requiring the
eigensystem of general matrices is uncommon, we decided not to provide such
facilities.
Consideration was also given to the inclusion of explicit facilities for LU
decomposition (as provided for example in the BLAS). LU decomposition is a
common first step for many operations. Thus making it available to the user
permits more rapid computation when several operations such as solving
equations and determinant evaluation are to be performed using
the same matrix. However, modern hardware is so fast that this would only
seem to be necessary for very large sets of equations and these are not the
target of these simple facilities. Moreover, adding explicit LU decomposition
introduces complexity for the user.
Consideration was also given to a fuller implementation of the BLAS. However,
this seems out of place in a language standard since it would be extremely long
and specialized. Such a fuller interface to the BLAS could be provided in
additional child packages. The goal here has been to provide convenient access
to simple type declarations and the more commonly required operations that are
not trivial for the user to program. These operations can of course be
implemented as a binding to an implementation of part of the BLAS.
!corrigendum G.3(01)
@dinsc
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.
!corrigendum G.3.1(01)
@dinsc
@i<@s8>
The generic library package Numerics.Generic_Real_Arrays has the following
declaration:
@xcode<@b
@b Real @b <@>;
package Ada.Numerics.Generic_Real_Arrays is
@b Pure(Generic_Real_Arrays);
-- @ft<@i>
@b Real_Vector @b (Integer @b <@>) @b Real'Base;
@b Real_Matrix @b (Integer @b <@>, Integer @b <@>) @b Real'Base;
-- @ft<@i>
-- @ft<@i>
@b "+" (Right : Real_Vector) @b Real_Vector;
@b "-" (Right : Real_Vector) @b Real_Vector;
@b "@b" (Right : Real_Vector) @b Real_Vector;
@b "+" (Left, Right : Real_Vector) @b Real_Vector;
@b "-" (Left, Right : Real_Vector) @b Real_Vector;
@b "*" (Left, Right : Real_Vector) @b Real'Base;
-- @ft<@i>
@b "*" (Left : Real'Base; Right : Real_Vector) @b Real_Vector;
@b "*" (Left : Real_Vector; Right : Real'Base) @b Real_Vector;
@b "/" (Left : Real_Vector; Right : Real'Base) @b Real_Vector;
-- @ft<@i>
@b Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) @b Real_Vector;
-- @ft<@i>
-- @ft<@i>
@b "+" (Right : Real_Matrix) @b Real_Matrix;
@b "-" (Right : Real_Matrix) @b Real_Matrix;
@b "@b" (Right : Real_Matrix) @b Real_Matrix;
@b Transpose (X : Real_Matrix) @b Real_Matrix;
@b "+" (Left, Right : Real_Matrix) @b Real_Matrix;
@b "-" (Left, Right : Real_Matrix) @b Real_Matrix;
@b "*" (Left, Right : Real_Matrix) @b Real_Matrix;
@b "*" (Left, Right : Real_Vector) @b Real_Matrix;
@b "*" (Left : Real_Vector; Right : Real_Matrix) @b Real_Vector;
@b "*" (Left : Real_Matrix; Right : Real_Vector) @b Real_Vector;
-- @ft<@i>
@b "*" (Left : Real'Base; Right : Real_Matrix) @b Real_Matrix;
@b "*" (Left : Real_Matrix; Right : Real'Base) @b Real_Matrix;
@b "/" (Left : Real_Matrix; Right : Real'Base) @b Real_Matrix;
-- @ft<@i>
@b Solve (A : Real_Matrix; X: Real_Vector) @b Real_Vector;
@b Solve (A, X : Real_Matrix) @b Real_Matrix;
@b Inverse (A : Real_Matrix) @b Real_Matrix;
@b Determinant (A : Real_Matrix) @b Real'Base;
-- @ft<@i>
@b Eigenvalues(A : Real_Matrix) @b Real_Vector;
@b Eigensystem(A : @b Real_Matrix;
Values : @b Real_Vector;
Vectors : @b Real_Matrix);
-- @ft<@i>
@b Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1)
@b Real_Matrix;
@b Ada.Numerics.Generic_Real_Arrays;>
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.
@xcode<@b "+" (Right : Real_Vector) @b Real_Vector;
@b "-" (Right : Real_Vector) @b Real_Vector;
@b "@b" (Right : Real_Vector) @b Real_Vector;>
@xindent
@xcode<@b "+" (Left, Right : Real_Vector) @b Real_Vector;
@b "-" (Left, Right : Real_Vector) @b Real_Vector;>
@xindent
@xcode<@b "*" (Left, Right : Real_Vector) @b Real'Base;>
@xindent
@xcode<@b "*" (Left : Real'Base; Right : Real_Vector) @b Real_Vector;>
@xindent
@xcode<@b "*" (Left : Real_Vector; Right : Real'Base) @b Real_Vector;
@b "/" (Left : Real_Vector; Right : Real'Base) @b Real_Vector;>
@xindent
@xcode<@b Unit_Vector (Index : Integer;
Order : Positive;
First : Integer := 1) @b Real_Vector;>
@xindent
First + Order - 1 or if First + Order - 1 @> Integer'Last.>
@xcode<@b "+" (Right : Real_Matrix) @b Real_Matrix;
@b "-" (Right : Real_Matrix) @b Real_Matrix;
@b "@b" (Right : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b Transpose (X : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b "+" (Left, Right : Real_Matrix) @b Real_Matrix;
@b "-" (Left, Right : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b "*" (Left, Right : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b "*" (Left, Right : Real_Vector) @b Real_Matrix;>
@xindent
@xcode<@b "*" (Left : Real_Vector; Right : Real_Matrix) @b Real_Vector;>
@xindent
@xcode<@b "*" (Left : Real_Matrix; Right : Real_Vector) @b Real_Vector;>
@xindent
@xcode<@b "*" (Left : Real'Base; Right : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b "*" (Left : Real_Matrix; Right : Real'Base) @b Real_Matrix;
@b "/" (Left : Real_Matrix; Right : Real'Base) @b Real_Matrix;>
@xindent
@xcode<@b Solve (A : Real_Matrix; X: Real_Vector) @b Real_Vector;>
@xindent
@xcode<@b Solve (A, X : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b Inverse (A : Real_Matrix) @b Real_Matrix;>
@xindent
@xcode<@b Determinant (A : Real_Matrix) @b Real'Base;>
@xindent
@xcode<@b Eigenvalues(A : Real_Matrix) @b Real_Vector;>
@xindent
@xcode<@b Eigensystem(A : @b Real_Matrix;
Values : @b Real_Vector;
Vectors : @b Real_Matrix);>
@xindent
@xcode<@b Unit_Matrix (Order : Positive;
First_1, First_2 : Integer := 1) @b Real_Matrix;>
@xindent Integer'Last or First_2 + Order - 1 @> Integer'Last.>
@i<@s8>
Accuracy requirements for the subprograms Solve, Inverse, Determinant,
Eigenvalues and Eigensystem 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).
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
@xindent
@i<@s8>
Implementations shall document any techniques used to reduce cancellation
errors such as extended precision arithmetic.
@i<@s8>
The nongeneric equivalent packages may, but need not, be actual
instantiations of the generic package for the appropriate predefined type.
@i<@s8>
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.
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.
The test that a matrix is symmetric may be performed by using the equality
operator to compare the relevant components.
!corrigendum G.3.2(01)
@dinsc
@i<@s8>
The generic library package Numerics.Generic_Complex_Arrays has the following
declaration:
@xcode<@b Ada.Numerics.Generic_Real_Arrays, Ada.Numerics.Generic_Complex_Types;
@b
@b Real_Arrays @b Ada.Numerics.Generic_Real_Arrays (<@>);
@b