Version 1.9 of ais/ai-00296.txt

Unformatted version of ais/ai-00296.txt version 1.9
Other versions for file ais/ai-00296.txt

!standard G (00)          03-07-23 AI95-00296/04
!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 including the use of LU decomposition. 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.
generic child packages for the determination of eigenvalues and eigenvectors are defined in G.3.3.
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; type Real_LU (<>) is private; -- LU decomposition of square matrix
-- 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 LU decomposition and related subprograms
function Decompose (A : Real_Matrix) return Real_LU; function Lower (LU: Real_LU) return Real_Matrix; function Upper (LU: Real_LU) return Real_Matrix; function Permutation (LU: Real_LU) return Real_Matrix;
-- Real_Matrix inversion and related operations
function Solve (A : Real_Matrix; X: Real_Vector) return Real_Vector; function Solve (LU : Real_LU; X: Real_Vector) return Real_Vector; function Solve (A, X : Real_Matrix) return Real_Matrix; function Solve (LU : Real_LU; X : Real_Matrix) return Real_Matrix;
function Inverse (A : Real_Matrix) return Real_Matrix; function Inverse (LU : Real_LU) return Real_Matrix; function Determinant (A : Real_Matrix) return Real'Base; function Determinant (LU : Real_LU) return Real'Base;
-- Other Real_Matrix operations
function Unit_Matrix (Order : Positive; First_1, First_2 : Integer := 1) return Real_Matrix;
private
-- implementation-defined
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.
Three 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 type Real_LU is a private type with unknown discriminants; values of this type represent real square matrices decomposed into lower and upper triangular matrices with possible row permutations.
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, 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 : 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. 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.
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 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, 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 : 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.
function Decompose (A : Real_Matrix) return Real_LU;
This function returns the LU decomposed form of the matrix A. Constraint_Error is raised if A'Length(1) and A'Length(2) are not equal.
function Lower (LU: Real_LU) return Real_Matrix;
This function returns the lower triangular matrix of the decomposed matrix LU. The diagonal components are 1.0. Both index ranges of the matrix result are A'Range(1) where A is the matrix of which LU is the decomposition.
function Upper (LU: Real_LU) return Real_Matrix;
This function returns the upper triangular matrix of the decomposed matrix LU. The index ranges of the matrix result are those of A where A is the matrix of which LU is the decomposition.
function Permutation (LU: Real_LU) return Real_Matrix;
This function returns a matrix whose components are zero except for one component in each row and column which is 1.0. It represents the row interchanges performed during the decomposition. Both index ranges of the result are A'Range(1) where A is the matrix of which LU is the decomposition.
function Solve (A : Real_Matrix; X: Real_Vector) return Real_Vector; function Solve (LU : Real_LU; X: Real_Vector) return Real_Vector;
These functions return a vector Y such that X is (nearly) equal to A * Y where A is either provided directly or is the matrix of which LU is the decomposition. 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; function Solve (LU : Real_LU; X : Real_Matrix) return Real_Matrix;
These functions return a matrix Y such that X is (nearly) equal to A * Y where A is either provided directly or is the matrix of which LU is the decomposition. 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; function Inverse (LU : Real_LU) return Real_Matrix;
These functions return a matrix B such that A * B is (nearly) the unit matrix where A is either provided directly or is the matrix of which LU is the decomposition. 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; function Determinant (LU : Real_LU) return Real'Base;
These functions return the determinant of the matrix A where A is either provided directly or is the matrix of which LU is the decomposition. Constraint_Error is raised if A'Length(1) is not equal to A'Length(2).
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. 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 Decompose, Lower, Upper, 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).
For operations involving an inner product, no requirements are specified in the relaxed mode. In the strict mode the accuracy 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 LU decomposition with row interchanges. This could be via a binding to an implementation of the so-called BLAS (Basic Linear Algebra Subprograms). Bindings to other facilities in the BLAS could be provided in child packages.
The implementation of the type Real_LU must retain information regarding the ranges of the original array since these are needed by Solve and Inverse for compatibility with the direct versions.
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.
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; type Complex_LU (<>) is private; -- LU decomposition of square matrix
-- 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 LU decomposition and related subprograms
function Decompose (A : Complex_Matrix) return Complex_LU; function Lower (LU: Complex_LU) return Complex_Matrix; function Upper (LU: Complex_LU) return Complex_Matrix; function Permutation (LU: Complex_LU) return Complex_Matrix;
-- Complex_Matrix inversion and related operations
function Solve (A : Complex_Matrix; X: Complex_Vector) return Complex_Vector; function Solve (LU : Complex_LU; X: Complex_Vector) return Complex_Vector; function Solve (A, X : Complex_Matrix) return Complex_Matrix; function Solve (LU : Complex_LU; X : Complex_Matrix) return Complex_Matrix;
function Inverse (A : Complex_Matrix) return Complex_Matrix; function Inverse (A : Complex_LU) return Complex_Matrix; function Determinant (A : Complex_Matrix) return Complex; function Determinant (A : Complex_LU) return Complex;
-- 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.
Three 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 type Complex_LU is a private type with unknown discriminants; values of this type represent complex square matrices decomposed into lower and upper triangular matrices with possible row permutations.
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.
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. The exception 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. The exception 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. 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.
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. The exception 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. The exception 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. The exception 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. The exception 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). The exception 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. The exception 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. 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).
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. 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.
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. 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, 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). 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_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). 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.
function Decompose (A : Complex_Matrix) return Complex_LU;
This function returns the LU decomposed form of the matrix A. Constraint_Error is raised if A'Length(1) and A'Length(2) are not equal.
function Lower (LU: Complex_LU) return Complex_Matrix;
This function returns the lower triangular matrix of the decomposed matrix LU. The diagonal components are 1.0. Both index ranges of the matrix result are A'Range(1) where A is the matrix of which LU is the decomposition.
function Upper (LU: Complex_LU) return Complex_Matrix;
This function returns the upper triangular matrix of the decomposed matrix LU. The index ranges of the matrix result are those of A where A is the matrix of which LU is the decomposition.
function Permutation (LU: Complex_LU) return Real_Matrix;
This function returns a real matrix whose components are zero except for one component in each row and column which is 1.0. It represents the row interchanges performed during the decomposition. Both index ranges of the result are A'Range(1) where A is the matrix of which LU is the decomposition.
function Solve (A : Complex_Matrix; X: Complex_Vector) return Complex_Vector; function Solve (LU : Complex_LU; X: Complex_Vector) return Complex_Vector;
These functions return a vector Y such that X is (nearly) equal to A * Y where A is either provided directly or is the matrix of which LU is the decomposition. 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; function Solve (LU : Complex_LU; X : Complex_Matrix) return Complex_Matrix;
These functions return a matrix Y such that X is (nearly) equal to A * Y where A is either provided directly or is the matrix of which LU is the decomposition. 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; function Inverse (LU : Complex_LU) return Complex_Matrix;
These functions return a matrix B such that A * B is (nearly) the unit matrix where A is either provided directly or is the matrix of which LU is the decomposition. 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; function Determinant (LU : Complex_LU) return Complex;
These functions return the determinant of the matrix A where A is either provided directly or is the matrix of which LU is the decomposition. Constraint_Error is raised if A'Length(1) is not equal to A'Length(2).
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). 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 Decompose, Lower, Upper, 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).
For operations involving an inner product, no requirements are specified in the relaxed mode. In the strict mode the accuracy 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 LU decomposition with row interchanges. This could be via a binding to an implementation of the so-called BLAS (Basic Linear Algebra Subprograms). Bindings to other facilities in the BLAS could be provided in child packages.
The implementation of the type Complex_LU must retain information regarding the ranges of the original array since these are needed by Solve and Inverse for compatibility with the direct versions.
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.
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:
generic package Ada.Numerics.Generic_Real_Arrays.Eigen is pragma Pure(Eigen);
function Values(A : Real_Matrix) return Real_Vector; function Values(LU : Real_LU) return Real_Vector;
procedure Values_And_Vectors(A : in Real_Matrix; Values : out Real_Vector; Vectors : out Real_Matrix); procedure Values_And_Vectors(LU : in Real_LU; Values : out Real_Vector; Vectors : out Real_Matrix);
end Ada.Numerics.Generic_Real_Arrays.Eigen;
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.
function Values(A : Real_Matrix) return Real_Vector; function Values(LU : Real_LU) return Real_Vector;
These functions return the eigenvalues of the symmetric matrix A as a vector sorted into order with the largest first. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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). If the matrix A is not symmetric then Argument_Error is raised.
procedure Values_And_Vectors(A : in Real_Matrix; Values : out Real_Vector; Vectors : out Real_Matrix); procedure Values_And_Vectors(LU : in Real_LU; Values : out Real_Vector; Vectors : out Real_Matrix);
These procedures compute both the eigenvalues and eigenvectors of the symmetric matrix A. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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. If the matrix A is not symmetric then Argument_Error is raised.
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:
generic package Ada.Numerics.Generic_Complex_Arrays.Eigen is pragma Pure(Eigen);
-- Eigensystem of a general real matrix
function Values(A : Real_Matrix) return Complex_Vector; function Values(LU : Real_LU) return Complex_Vector;
procedure Values_And_Vectors(A : in Real_Matrix; Values : out Complex_Vector; Vectors : out Complex_Matrix); procedure Values_And_Vectors(LU : in Real_LU; Values : out Complex_Vector; Vectors : out Complex_Matrix);
-- Eigensystem of a Hermitian matrix
function Values(A : Complex_Matrix) return Real_Vector; function Values(LU : Complex_LU) return Real_Vector;
procedure Values_And_Vectors(A : in Complex_Matrix; Values : out Real_Vector; Vectors : out Complex_Matrix); procedure Values_And_Vectors(LU : in Complex_LU; Values : out Real_Vector; Vectors : out Complex_Matrix);
-- Eigensystem of a general complex matrix
function Values(A : Complex_Matrix) return Complex_Vector; function Values(LU : Complex_LU) return Complex_Vector;
procedure Values_And_Vectors(A : in Complex_Matrix; Values : out Complex_Vector; Vectors : out Complex_Matrix); procedure Values_And_Vectors(LU : in Complex_LU; Values : out Complex_Vector; Vectors : out Complex_Matrix);
end Ada.Numerics.Generic_Complex_Arrays.Eigen;
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.
The effect of the subprograms is as described below.
function Values(A : Real_Matrix) return Complex_Vector; function Values(LU : Real_LU) return Complex_Vector;
These functions return the eigenvalues of the matrix A as a vector sorted by order of modulus with the largest first. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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).
procedure Values_And_Vectors(A : in Real_Matrix; Values : out Complex_Vector; Vectors : out Complex_Matrix); procedure Values_And_Vectors(LU : in Real_LU; Values : out Complex_Vector; Vectors : out Complex_Matrix);
These procedures compute both the eigenvalues and eigenvectors of the matrix A. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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.
function Values(A : Complex_Matrix) return Real_Vector; function Values(LU : Complex_LU) return Real_Vector;
These functions return the eigenvalues of the Hermitian matrix A as a vector sorted into order with the largest first. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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 Argument_Error is raised if the matrix is not Hermitian.
procedure Values_And_Vectors(A : in Complex_Matrix; Values : out Real_Vector; Vectors : out Complex_Matrix); procedure Values_And_Vectors(LU : in Complex_LU; Values : out Real_Vector; Vectors : out Complex_Matrix);
These procedures compute both the eigenvalues and eigenvectors of the Hermitian matrix A. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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 Argument_Error is raised if the matrix is not Hermitian.
function Values(A : Complex_Matrix) return Complex_Vector; function Values(LU : Complex_LU) return Complex_Vector;
These functions return the eigenvalues of the matrix A as a vector sorted by order of modulus with the largest first. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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).
procedure Values_And_Vectors(A : in Complex_Matrix; Values : out Complex_Vector;
Vectors : out Complex_Matrix); procedure Values_And_Vectors(LU : in Complex_LU; Values : out Complex_Vector; Vectors : out Complex_Matrix);
These procedures compute both the eigenvalues and eigenvectors of the matrix A. The matrix A is either provided directly or is the matrix of which LU is the decomposition. 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.
Implementation Advice
The test that a matrix is symmetric may be performed by using the equality operator to compare the relevant components. The test that a matrix is Hermitian may similarly use the equality operator to compare the real components and negation followed by equality to compare the imaginary components (see G.2.1).
?? I hope that is good enough. I am not sure whether the RM implies that negation is exact. If it doesn't then the paragraph above needs changing or G.2.1 needs changing.
AARM Note
In this child package, Generic_Complex_Arrays.Eigen, the full view of the type Real_LU is not available. However, the details may be obtained by calling the functions Lower, Upper and Permutation in the formal package Real_Arrays. This is the real reason for declaring these subprograms although they might be helpful to the inquisitive user.
!example
This example illustrates the use of the LU decomposed form to avoid unnecessary computation.
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; with Ada.Numerics.Real_Arrays.Eigen; procedure An_Example is N : Integer := ... Mat : Real_Array(1 .. N, 1 .. N); P, Q, R : Real_Vector(1 .. N); X, Y, Z : Real_Vector(1 .. N); -- compute some value for Mat and values for vectors P, Q and R;
begin LUMat : constant Real_LU := Decompose(Mat); -- Now compute various properties of Mat D : Float; Inverse_Mat : Real_Array(1 .. N); Mat_Eigenvalues : Real_Array(1 .. N); ... D := Determinant(LUMat); Inverse_Mat := Invert(LUMat); X := Solve(LUMat, P); Y := Solve(LUMat, Q); Z := Solve(LUMat, R); Mat_Eigenvalues := Eigen.Values(LUMat); ... end; end An_Example;
!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 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.
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 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.
In addition a private type has been added to hold the LU decomposed form of a matrix. This enables repeated calls of Solve etc to be done efficiently by avoiding unnecessary repeated decomposition. Further overloadings of the various subprograms using the LU form are thus provided. This private type has unknown discriminants and thus objects of this type must be initialized with a call of Decompose. The use of a private type prevents the user from making errors by confusing the original and LU forms of a matrix.
The private type has unknown discriminants. The full type must retain the information regarding the bounds of the original array since it is needed for some subprograms. The functions Lower, Upper and Permutation give access to the details of the decomposition. Note that the permutation describing the row interchanges is given in the rather unorthodox form of a real matrix. It would be more conventional just to provide an integer vector defining the permutation but there is no convenient type available and declaring an integer vector type just for this seems tiresome. Moreover, the matrix form of the permutation enables the original array to be reconstructed simply by multiplication.
?? I hope this is true. Someone better check this carefully.
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 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 subprograms 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 functions for the eigenvectors).
In the case of eigenvectors of matrices which are neither symmetric nor Hermitian no requirement has been placed on the independence or modulus of the vectors.
All the eigenvalue and vector subprograms also have overloadings for using the LU decomposed forms. A curiosity arises in the case of the eigenvalues of a general real matrix because this needs a parameter of type Real_LU. However, the full view of this type is not available in the complex eigenvalue package and so the decomposed form has to be accessed via the subprograms Lower, Upper and Permutation.
Consideration was 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 easy for the user to program. These operations can of course be implemented as a binding to an implementation of part of the BLAS.
!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.

****************************************************************

From: John Barnes
Sent: Thursday, July 24, 2003  2:09 AM

Version /03 is as /02 dated 03-05-22 but has vector product removed
and other minor changes agreed at the ARG meeting in Toulouse. It does not
include LU decomposition.

****************************************************************

From: John Barnes
Sent: Thursday, July 24, 2003  2:09 AM

Version /04 is as /03 dated 03-07-22 but includes LU decomposition.

****************************************************************

From: John Barnes
Sent: Friday, July 25, 2003  2:29 AM

I have created two new versions of this AI and Randy has put
them on the database (he tells me I had some spurious CRs
which he had to edit and so he warns that perhaps they are
not perfect in that respect).

The first (version 03) incorporates various minor changes
agreed at Toulouse but does not include LU decomposition.
The second (version 04) includes LU decomposition.

The ZIP file only contains 04. Randy tells me so you will
have to go to the database directly for 03. (I always use
the database directly so maybe that's no problem.)


You may recall (or then you may not) that LU decomposition
enables a square matrix to be held in a form which is faster
for operations such as invert, solve and eigenvalues.  It is
of value when several operations have to be performed on the
same matrix because time is saved by not repeating some
common operations.

This was important when hardware was slow. But hardware is
so fast now that it can hardly matter in most circumstances.
Perhaps it might matter in a tight loop when the same matrix
is being used to solve a set of more than 100 equations
repeatedly. But does this ever happen?  Seems doubtful and
if it did a real numerical analyst needs to be called in
anyway.

Note that large banded matrices do occur in differential
equations but that is specialized stuff anyway which we are
not addressing.

Moreover adding LU adds complexity to the RM and an
irritating visibility problem explained in the AI
discussion.

Since our goal here I believe is to provide some simple
useful facilities for the common programmer (while still
providing a base which can be used through child packages to
provide other stuff) it all now seems too much to me.

As a programmer nearly half a century ago, I used to invert
lots of smallish matrices and never used LU decomposition
(maybe it hadn't been discovered/invented).  Maybe there is
some partial analogy with all those wonderful algorithms for
dealing with large matrices held on several magnetic tapes.
All in the past surely. Of course LU has some nice
mathematics behind it but not every Ada programmer needs to
be exposed to such things.

I am sure that any professional numerical analyst will say
that LU is a must have. But that could be more out of
professional pride, habit and rectitude rather than
practical need.

I am also sure that if presented with the sort of numerical
problems I dabble with from time to time I would never
bother with the LU forms because the straightforward ones
would do the job adequately. Another point is that the
existence of the LU forms might make life harder for
software maintenance by obscuring what is going on.

Perhaps I should add that my colleagues in the UK who have
looked at this are moving to the view that LU is
unnecessary. And discussions with Pascal before he went on
vacation were going that way as well.

We can discuss it in Sydney. But any thoughts before then
might be welcome. I suppose I will have to do the wording
for downward closures now - groan.

****************************************************************


Questions? Ask the ACAA Technical Agent