CVS difference for ais/ai-00296.txt

Differences between 1.3 and version 1.4
Log of other versions for file ais/ai-00296.txt

--- ais/ai-00296.txt	2003/01/24 04:14:27	1.3
+++ ais/ai-00296.txt	2003/06/05 00:34:14	1.4
@@ -1,4 +1,4 @@
-!standard G (00)                                        03-02-21  AI95-00296/01
+!standard G (00)                                        03-05-22  AI95-00296/02
 !class amendment 02-06-07
 !status received 02-06-07
 !priority Medium
@@ -7,8 +7,8 @@
 
 !summary
 
-The vector and matrix operations in ISO/IEC 13813 are added to Ada.Numerics in
-Annex G.
+The vector and matrix operations in ISO/IEC 13813 plus related operations
+are added to Ada.Numerics in Annex G.
 
 !problem
 
@@ -33,11 +33,9 @@
 They are included as a new subclause G.3 in order to avoid excessive
 renumbering of other clauses.
 
-However, if rearrangement of the Annex were deemed acceptable then it might be
-better to put the real arrays package in G.1 (any other non-complex numeric
-library packages might then also be in G.1), to renumber G.1 (Complex
-Arithmetic) as G.2 and G.2 (Numeric Performance Requirements) as G.3. The
-complex arrays package could then be G.2.5.
+In addition to the main facilities of 13813, these packages also include
+operations for matrix inversion and determinants. Furthermore, child packages
+provide subprograms for the determination of eigenvalues and eigenvectors.
 
 !wording
 
@@ -49,9 +47,12 @@
 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 predefind package Numerics (see A.5).
+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.
 
-??? do we need non-generic equivalents ???
+Generic child packages for the determination of eigenvalues and eigenvectors
+are defined in G.3.3.
 
 G.3.1 Real Vectors and Matrices
 
@@ -62,36 +63,33 @@
 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
+   -- Real_Vector additive 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_Vector;
-   function "/"   (Left, Right : Real_Vector) return Real_Vector;
-   function "**"  (Left        : Real_Vector;
-                   Right       : Integer)     return Real_Vector;
 
-   function "*"   (Left, Right : Real_Vector) return Real'Base;
+   function "+"  (Left, Right : Real_Vector) return Real_Vector;
+   function "-"  (Left, Right : Real_Vector) return Real_Vector;
 
    -- 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;
+   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;
+
+   -- Inner, Outer and Vector products
+
+   function "*" (Left, Right : Real_Vector) return Real'Base;
+   function "*" (Left, Right : Real_Vector) return Real_Matrix;
+   function "*" (Left, Right : Real_Vector) return Real_Vector;
 
    -- Other Real_Vector operations
 
@@ -99,8 +97,6 @@
                          Order : Positive;
                          First : Integer := 1) return Real_Vector;
 
-   -- Subprograms for Real_Matrix Types
-
    -- Real_Matrix arithmetic operations
 
    function "+"       (Right : Real_Matrix) return Real_Matrix;
@@ -111,25 +107,37 @@
    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;
+   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;
+   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;
 
    -- Other Real_Matrix operations
 
-   function Identity_Matrix (Order            : Positive;
-                             First_1, First_2 : Integer := 1) return Real_Matrix;
+   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
@@ -138,63 +146,66 @@
 a matrix with components of type Real; it is defined as an unconstrained,
 two-dimensional array with indices of type Integer.
 
-?? maybe indexes ??
-
 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. These are the standard mathematical
-operations for vector identity, negation and absolute value respectively. The
-bounds of the result are those of Right.
+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;
-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.
-These are the standard mathematical operations for vector addition,
-subtraction, multiplication and division respectively. The bounds of the result
-are those of the left operand. The exception Constraint_Error is raised if
-Left'Length is not equal to Right'Length.
+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 : Integer) return Real_Vector;
+function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector;
 
-This operation returns the result of applying the corresponding operation "**"
-of the type Real with integer power Right to each component of Left. The bounds
-of the result are those of the left operand.
+This operation returns the result of multiplying each component of Right by the
+scalar Left using the "*" operation of the type Real. The index range of the
+result is Right'Range.
+
+function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector;
+function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector;
+
+Each operation returns the result of applying the corresponding operation of
+the type Real to each component of Left and to the scalar Right. The index
+range of the result is Left'Range.
 
-function "*" (Left, Right  : Real_Vector) return Real'Base;
+function "*" (Left, Right : Real_Vector) return Real'Base;
 
-This operation returns the inner (dot) product of Left and Right. The exception
+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; an accuracy requirement is not specified.
-Constraint_Error is raised if an intermediate result is outside the range of
-Real'Base even though the mathematical final result would not be.
+operation involves an inner product.
 
-function "*" (Left  : Real'Base; Right : Real_Vector) return Real_Vector;
+function "*" (Left, Right : Real_Vector) 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. This is the standard
-mathematical operation for scaling a vector Right by a real number Left. The
-bounds of the result are those of the right operand.
+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'Base) return Real_Vector;
-function "/" (Left  : Real_Vector; Right : Real'Base) 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 to the scalar Right. These are
-standard mathematical operations for scaling a vector Left by a real number
-Right. The bounds of the result are those of the left operand.
+This operation returns the vector product of Left and Right. The exception
+Constraint_Error is raised if Left'Length and Right'Length are not equal to 3.
+The index range of the result is Left'Range.
+
+?? I am not sure about vector product. It only exists in three dimensions and
+is really a skew-symmetric tensor.
 
 function Unit_Vector (Index : Integer;
                       Order : Positive;
@@ -203,97 +214,141 @@
 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. This function is
-exact.
+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. These are the standard mathematical
-operations for matrix identity, negation and absolute value respectively. The
-bounds of the result are those of Right.
+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 index ranges of the
-result are X'Range(2) and X'Range(1) (first and second index respectively).
-This function is exact.
+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. These
-are also the standard mathematical operations for matrix addition and
-subtraction. The bounds of the result are those of Left. The exception
+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 index ranges of the result are Left'Range(1) and
-Right'Range(2) (first and second index respectively). The exception
+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 an inner product; an accuracy requirement is not
-specified. Constraint_Error is raised if an intermediate result is outside the
-range of Real'Base even though the mathematical final result would not be.
-
-function "*" (Left, Right : Real_Vector) return Real_Matrix;
+This operation involves inner products.
 
-This operation provides the standard mathematical operation for multiplication
-of a (column) vector Left by a (row) vector Right. The index ranges of the
-matrix result are Left'Range and Right'Range (first and second index
-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 an inner
-product; an accuracy requirement is not specified. Constraint_Error is raised
-if an intermediate result is outside the range of Real'Base even though the
-mathematical final result would not be.
+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 an inner
-product; an accuracy requirement is not specified. Constraint_Error is raised
-if an intermediate result is outside the range of Real'Base even though the
-mathematical final result would not be.
+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. This is the standard
-mathematical operation for scaling a matrix Right by a real number Left. The
-index ranges of the matrix result are those of Right.
+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;
+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. These are
-standard mathematical operations for scaling a matrix Left by a real number
-Right. The index ranges of the matrix result are those of Left.
+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) 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 Identity_Matrix (Order            : Positive;
+function Unit_Matrix (Order            : Positive;
                           First_1, First_2 : Integer := 1) return Real_Matrix;
 
-This function returns a square "identity matrix" with Order**2 components and
+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. This
-function is exact.
+First_1 + Order - 1 > Integer'Last or First_2 + Order - 1 > Integer'Last.
 
+Implementation Requirements
 
+Accuracy requirements for the functions Solve, Inverse and Determinant are
+implementation defined.
+
+For operations not involving an inner product, the accuracy
+requirements are those of the corresponding operations of the type Real in
+both the strict mode and the relaxed mode (see G.2).
+
+?? This doesn't really cover vector product but I might delete it anyway.
+
+For operations involving an inner product, no requirements are specified in
+the relaxed mode. In the strict mode the acuracy should be at least that of
+the canonical implementation of multiplication and addition using the
+corresponding operations of type Real'Base and performing the cumulative
+addition using ascending indices. Implementations shall document any
+techniques used to reduce cancellation errors such as extended precision
+arithmetic.
+
+Implementation Permissions
+
+The nongeneric equivalent packages may, but need not, be actual
+instantiations of the generic package for the appropriate predefined type.
+
+Implementation Advice
+
+Implementations should implement the Solve and Inverse functions using
+established techniques such as Gauss-Seidel with Interchanges.
+
+It is not the intention that any special provision should be made to
+determine whether a matrix is ill-conditioned or not. The naturally occuring
+division by zero or overflow which will result from executing these
+functions with an ill-conditioned matrix and thus raise Constraint_Error is
+sufficient.
+
+
 G.3.2 Complex Vectors and Matrices
 
 The generic library package Numerics.Generic_Complex_Arrays has the following
@@ -315,8 +370,6 @@
    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;
@@ -331,18 +384,19 @@
    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 "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;
+                                                    return Complex_Vector;
    function Compose_From_Polar (Modulus, Argument : Real_Vector;
                                 Cycle             : Real'Base)
-      return Complex_Vector;
+                                                    return Complex_Vector;
 
-   -- Complex_Vector arithmetic operations
+   -- Complex_Vector additive operations
 
    function "+"       (Right  : Complex_Vector) return Complex_Vector;
    function "-"       (Right  : Complex_Vector) return Complex_Vector;
@@ -350,14 +404,9 @@
 
    function "+"  (Left, Right : 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_Vector;
-   function "**" (Left  : Complex_Vector;
-                  Right : Integer)              return Complex_Vector;
 
-   function "*"  (Left, Right : Complex_Vector) return Complex;
 
-   -- Mixed Real_Vector and Complex_Vector arithmetic operations
+   -- Mixed Real_Vector and Complex_Vector additive operations
 
    function "+" (Left  : Real_Vector;
                  Right : Complex_Vector) return Complex_Vector;
@@ -367,17 +416,6 @@
                  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_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
 
@@ -395,14 +433,30 @@
    function "/" (Left  : Complex_Vector;
                  Right : Real'Base)      return Complex_Vector;
 
+   -- Inner, Outer and Vector products
+
+   function "*" (Left, Right : Complex_Vector) return Complex;
+   function "*" (Left  : Real_Vector;    Right : Complex_Vector) return Complex;
+   function "*" (Left  : Complex_Vector; Right : Real_Vector)    return Complex;
+
+   function "*" (Left, Right : Complex_Vector) return Complex_Matrix;
+   function "*" (Left  : Real_Vector;
+                 Right : Complex_Vector) return Complex_Matrix;
+   function "*" (Left  : Complex_Vector;
+                 Right : Real_Vector)    return Complex_Matrix;
+
+   function "*" (Left, Right : Complex_Vector) return Complex_Vector;
+   function "*" (Left  : Real_Vector;
+                 Right : Complex_Vector) return Complex_Vector;
+   function "*" (Left  : Complex_Vector;
+                 Right : Real_Vector)    return Complex_Vector;
+
    -- 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;
@@ -417,17 +471,18 @@
    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 "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;
+                                                    return Complex_Matrix;
    function Compose_From_Polar (Modulus, Argument : Real_Matrix;
                                 Cycle             : Real'Base)
-      return Complex_Matrix;
+                                                    return Complex_Matrix;
 
    -- Complex_Matrix arithmetic operations
 
@@ -440,8 +495,6 @@
    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;
@@ -463,10 +516,6 @@
                  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;
@@ -491,23 +540,791 @@
    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;
+
+??  I have not added mixed real and complex for Solve.
+
    -- Other Complex_Matrix operations
 
-   function Identity_Matrix (Order            : Positive;
-                             First_1, First_2 : Integer := 1)
-      return Complex_Matrix;
+   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 and Imaginary types exported by Numerics.Complex_Types are
+systematically substituted for Complex and Imaginary, throughout. Nongeneric
+equivalents for each of the other predefined floating point types are
+defined similarly, with the names Numerics.Short_Complex_Arrays,
+Numerics.Long_Complex_Arrays, etc.
+
+?? I based this text on that in G.1.2(9/1). However this clause (ie G.3.2)
+doesn't seem to use Imaginary but I left it in.
+
+Two types are defined and exported by Ada.Numerics.Generic_Complex_Arrays.
+The composite type Complex_Vector is provided to represent a vector with
+components of type Complex; it is defined as an unconstrained
+one-dimensional array with an index of type Integer. The composite type
+Complex_Matrix is provided to represent a matrix with components of type
+Complex; it is defined as an unconstrained, two-dimensional array with
+indices of type Integer.
+
+The effect of the various subprograms is as described below. In many cases
+they are described in terms of corresponding scalar operations in
+Numerics.Generic_Complex_Types. Any exception raised by those operations is
+propagated by the array subprogram. Moreover any constraints on the parameters
+and the accuracy of the result for each individual component is as defined for
+the scalar operation.
+
+In the case of those operations which are defined to involve an inner product,
+Constraint_Error may be raised if an intermediate result has a component
+outside the range of Real'Base even though the final mathematical result
+would not.
+
+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 and if X'Length is not equal to 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;
+
+function "+" (Left  : Real_Vector;
+              Right : Complex_Vector) return Complex_Vector;
+function "+" (Left  : Complex_Vector;
+              Right : Real_Vector)    return Complex_Vector;
+function "-" (Left  : Real_Vector;
+              Right : Complex_Vector) return Complex_Vector;
+function "-" (Left  : Complex_Vector;
+              Right : Real_Vector)    return Complex_Vector;
+
+Each operation returns the result of applying the corresponding operation
+in Numerics.Generic_Complex_Types to each component of Left and the
+matching component of Right. The index range of the result is Left'Range.
+The exception Constraint_Error is raised if Left'Length is not equal to
+Right'Length.
+
+function "*" (Left : Complex; Right : Complex_Vector) return Complex_Vector;
+
+This operation returns the result of multiplying each component of Right by
+the complex number Left using the appropriate operation "*" in
+Numerics.Generic_Complex_Types. The index range of the result is Right'Range.
+
+function "*" (Left : Complex_Vector; Right : Complex) return Complex_Vector;
+function "/" (Left : Complex_Vector; Right : Complex) return Complex_Vector;
+
+Each operation returns the result of applying the corresponding operation in
+Numerics.Generic_Complex_Types to each component of the vector Left and the
+complex number Right. The index range of the result is Left'Range.
+
+function "*" (Left : Real'Base; Right : Complex_Vector) return Complex_Vector;
+
+This operation returns the result of multiplying each component of Right by
+the real number Left using the appropriate operation "*" in
+Numerics.Generic_Complex_Types. The index range of the result is Right'Range.
+
+function "*" (Left : Complex_Vector; Right : Real'Base) return Complex_Vector;
+function "/" (Left : Complex_Vector; Right : Real'Base) return Complex_Vector;
+
+Each operation returns the result of applying the corresponding operation in
+Numerics.Generic_Complex_Types to each component of the vector Left and the
+real number Right. The index range of the result is Left'Range.
+
+function "*" (Left, Right : Complex_Vector) return Complex;
+function "*" (Left : Real_Vector;    Right : Complex_Vector) return Complex;
+function "*" (Left : Complex_Vector; Right : Real_Vector)    return Complex;
+
+Each operation returns the inner product of Left and Right. The exception
+Constraint_Error is raised if Left'Length is not equal to Right'Length.
+These operations involve an inner product.
+
+function "*" (Left, Right : Complex_Vector) return Complex_Matrix;
+function "*" (Left  : Real_Vector;
+              Right : Complex_Vector) return Complex_Matrix;
+function "*" (Left  : Complex_Vector;
+              Right : Real_Vector)    return Complex_Matrix;
+
+Each operation returns the outer product of a (column) vector Left by a (row)
+vector Right using the appropriate operation "*" in
+Numerics.Generic_Complex_Types for computing the individual components.
+The first and second index ranges of the matrix result are Left'Range and
+Right'Range respectively.
+
+function "*" (Left, Right : Complex_Vector) return Complex_Vector;
+function "*" (Left  : Real_Vector;
+              Right : Complex_Vector) return Complex_Vector;
+function "*" (Left  : Complex_Vector;
+              Right : Real_Vector)    return Complex_Vector;
+
+Each operation returns the vector product of Left and Right. The exception
+Constraint_Error is raised if Left'Length and Right'Length are not equal to 3.
+The index range of the result is Left'Range.
+
+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
+X'Length(2) is not equal to Re'Length(2) and if X'Length(1) is not equal to
+Im'Length(1) or X'Length(2) is not equal to Im'Length(2).
+
+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  : Complex_Vector;
+              Right : Complex_Matrix) return Complex_Vector;
+
+This operation provides the standard mathematical operation for multiplication
+of a (row) vector Left by a matrix Right. The index range of the (row) vector
+result is Right'Range(2). The exception Constraint_Error is raised if
+Left'Length is not equal to Right'Length(1). This operation involves inner
+products.
+
+function "*" (Left  : Complex_Matrix;
+              Right : Complex_Vector) return Complex_Vector;
+
+This operation provides the standard mathematical operation for multiplication
+of a matrix Left by a (column) vector Right. The index range of the (column)
+vector result is Left'Range(1). The exception Constraint_Error is raised if
+Left'Length(2) is not equal to Right'Length. This operation involves inner
+products.
+
+function "+" (Left  : Real_Matrix;
+              Right : Complex_Matrix) return Complex_Matrix;
+function "+" (Left  : Complex_Matrix;
+              Right : Real_Matrix)    return Complex_Matrix;
+function "-" (Left  : Real_Matrix;
+              Right : Complex_Matrix) return Complex_Matrix;
+function "-" (Left  : Complex_Matrix;
+              Right : Real_Matrix)    return Complex_Matrix;
+
+Each operation returns the result of applying the corresponding operation in
+Numerics.Generic_Complex_Types to each component of Left and the matching
+component of Right. The index ranges of the result are those of Left. The
+exception Constraint_Error is raised if Left'Length(1) is not equal to
+Right'Length(1) or Left'Length(2) is not equal to Right'Length(2).
+
+function "*" (Left  : Real_Matrix;
+              Right : Complex_Matrix) return Complex_Matrix;
+function "*" (Left  : Complex_Matrix;
+              Right : Real_Matrix)    return Complex_Matrix;
+
+Each operation provides the standard mathematical operation for matrix
+multiplication. The first and second index ranges of the result are
+Left'Range(1) and Right'Range(2) respectively. The exception
+Constraint_Error is raised if Left'Length(2) is not equal to Right'Length(1).
+These operations involve inner products.
+
+function "*" (Left  : Real_Vector;
+              Right : Complex_Matrix) return Complex_Vector;
+function "*" (Left  : Complex_Vector;
+              Right : Real_Matrix)    return Complex_Vector;
+
+Each operation provides the standard mathematical operation for multiplication
+of a (row) vector Left by a matrix Right. The index range of the (row) vector
+result is Right'Range(2). The exception Constraint_Error is raised if
+Left'Length is not equal to Right'Length(1). These operations involve inner
+products.
+
+function "*" (Left  : Real_Matrix;
+              Right : Complex_Vector) return Complex_Vector;
+function "*" (Left  : Complex_Matrix;
+              Right : Real_Vector)    return Complex_Vector;
+
+Each operation provides the standard mathematical operation for multiplication
+of a matrix Left by a (column) vector Right. The index range of the (column)
+vector result is Left'Range(1). The exception Constraint_Error is raised if
+Left'Length(2) is not equal to Right'Length. These operations involve inner
+products.
+
+function "*" (Left : Complex; Right : Complex_Matrix) return Complex_Matrix;
+
+This operation returns the result of multiplying each component of Right by
+the complex number Left using the appropriate operation "*" in
+Numerics.Generic_Complex_Types. The index ranges of the result are those of
+Right.
+
+function "*" (Left : Complex_Matrix; Right : Complex) return Complex_Matrix;
+function "/" (Left : Complex_Matrix; Right : Complex) return Complex_Matrix;
+
+Each operation returns the result of applying the corresponding operation in
+Numerics.Generic_Complex_Types to each component of the matrix Left and the
+complex number Right. The index ranges of the result are those of Left.
+
+function "*" (Left : Real'Base; Right : Complex_Matrix) return Complex_Matrix;
+
+This operation returns the result of multiplying each component of Right by
+the real number Left using the appropriate operation "*" in
+Numerics.Generic_Complex_Types. The index ranges of the result are those of
+Right.
+
+function "*" (Left : Complex_Matrix; Right : Real'Base) return Complex_Matrix;
+function "/" (Left : Complex_Matrix; Right : Real'Base) return Complex_Matrix;
+
+Each operation returns the result of applying the corresponding operation in
+Numerics.Generic_Complex_Types to each component of the matrix Left and the
+scalar Right. The index ranges of the result are those of Left.
+
+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) 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 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 Solve, Inverse and Determinant are
+implementation defined.
+
+For operations not involving an inner product, the accuracy requirements are
+those of the corresponding operations of the type Real'Base and Complex in both
+the strict mode and the relaxed mode (see G.2).
+
+?? This doesn't really cover vector product but I might delete it anyway.
+
+For operations involving an inner product, no requirements are specified in
+the relaxed mode. In the strict mode the acuracy should be at least that of
+the canonical implementation of multiplication and addition using the
+corresponding operations of type Real'Base and Complex and performing the
+cumulative addition using ascending indices. Implementations shall document
+any techniques used to reduce cancellation errors such as extended precision
+arithmetic.
+
+Implementation Permissions
+
+The nongeneric equivalent packages may, but need not, be actual
+instantiations of the generic package for the appropriate predefined type.
+
+Although many operations are defined in terms of operations from
+Numerics.Generic_Complex_Types, they need not be implemented by calling those
+operations provided that the effect is the same.
+
+Implementation Advice
+
+Implementations should implement the Solve and Inverse functions using
+established techniques such as Gauss-Seidel with Interchanges.
+
+It is not the intention that any special provision should be made to
+determine whether a matrix is ill-conditioned or not. The naturally occuring
+division by zero or overflow which will result from executing these
+functions with an ill-conditioned matrix and thus raise Constraint_Error is
+sufficient.
+
+Implementations should not perform operations on mixed complex and real operands
+by first converting the real operand to complex. See G.1.1(56,57).
+
+G.3.3 Eigensystems of Matrices
+
+Eigenvalues and eigenvectors of a real symmetric matrix are computed using
+subprograms in the child package Ada.Numerics.Generic_Real_Arrays.Eigen
+whose declaration is:
+
+generic
+package Ada.Numerics.Generic_Real_Arrays.Eigen is
+   pragma Pure(Eigen);
+
+   function Values(A : Real_Matrix) return Real_Vector;
+
+   procedure Values_And_Vectors(A       : in  Real_Matrix;
+                                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;
+
+This function returns the eigenvalues of the symmetric matrix A as a vector
+sorted into order with the largest first. The exception Constraint_Error is
+raised if A'Length(1) is not equal to A'Length(2). The index range of the
+result is A'Range(1).
+
+procedure Values_And_Vectors(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 Values. The out parameter Vectors is a matrix whose columns are the
+eigenvectors of the matrix A. The order of the columns corresponds to the
+order of the eigenvalues. The eigenvectors are normalized and mutually
+orthogonal (they are orthonormal),including when there are repeated
+eigenvalues. The exception Constraint_Error is raised if A'Length(1) is not
+equal to A'Length(2). The index ranges of the parameter Vectors are those of A.
+
+For both subprograms, if the matrix A is not symmetric then Constraint_Error is
+raised.
+
+?? Or we could define the matrix just in terms of the upper diagonal or we
+could make it symmetric by taking the arithmetic mean of the matrix and its
+transpose.
+
+Eigenvalues and eigenvectors of a general real or complex matrix are computed
+using subprograms in the child package
+Ada.Numerics.Generic_Complex_Arrays.Eigen whose declaration is:
+
+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;
+
+   procedure Values_And_Vectors(A       : in  Real_Matrix;
+                                Values  :  out Complex_Vector;
+                                Vectors :  out Complex_Matrix);
+
+   -- Eigensystem of a Hermitian matrix
+
+   function Values(A : Complex_Matrix) return Real_Vector;
+
+   procedure Values_And_Vectors(A       : in  Complex_Matrix;
+                                Values  :  out Real_Vector;
+                                Vectors :  out Complex_Matrix);
+
+   -- Eigensystem of a general complex matrix
+
+   function Values(A : Complex_Matrix) return Complex_Vector;
+
+   procedure Values_And_Vectors(A       : in  Complex_Matrix;
+                                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.
+
+?? I didn't mention Real'Base, Complex etc because they are not used in the
+spec.  Should I??
+
+The effect of the subprograms is as described below.
+
+function Values(A : Real_Matrix) return Complex_Vector;
+
+This function returns the eigenvalues of the matrix A as a vector sorted by
+order of modulus with the largest first. The exception Constraint_Error is
+raised if A'Length(1) is not equal to A'Length(2). The index range of the
+result is A'Range(1).
+
+procedure Values_And_Vectors(A       : in  Real_Matrix;
+                             Values  :  out Complex_Vector;
+                             Vectors :  out Complex_Matrix);
+
+This procedure computes both the eigenvalues and eigenvectors of the
+matrix A. The out parameter Values is the same as that obtained by calling the
+function Values. The out parameter Vectors is a matrix whose columns are the
+eigenvectors of the matrix A. The order of the columns corresponds to the
+order of the eigenvalues. The exception Constraint_Error is
+raised if A'Length(1) is not equal to A'Length(2). The index ranges of the
+parameter Vectors are those of A.
+
+function Values(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. The exception Constraint_Error is
+raised if A'Length(1) is not equal to A'Length(2). The index range of the
+result is A'Range(1). The exception Constraint_Error is raised if the matrix
+is not hermitian.
+
+procedure Values_And_Vectors(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 Values. The out parameter Vectors is a matrix whose columns are the
+eigenvectors of the matrix A. The order of the columns corresponds to the
+order of the eigenvalues. The eigenvectors are mutually orthonormal,
+including when there are repeated eigenvalues. The exception Constraint_Error
+is raised if A'Length(1) is not equal to A'Length(2). The index ranges of the
+parameter Vectors are those of A. The exception Constraint_Error is raised if
+the matrix is not hermitian.
+
+?? Or we could define the matrix just in terms of the upper diagonal or we
+could make it hermitian by taking the arithmetic mean of the matrix and its
+conjugate transpose.
+
+function Values(A : Complex_Matrix) return Complex_Vector;
+
+This function returns the eigenvalues of the matrix A as a vector sorted by
+order of modulus with the largest first. The exception Constraint_Error is
+raised if A'Length(1) is not equal to A'Length(2). The index range of the
+result is A'Range(1).
+
+procedure Values_And_Vectors(A       : in  Complex_Matrix;
+                             Values  :  out Complex_Vector;
+                             Vectors :  out Complex_Matrix);
+
+This procedure computes both the eigenvalues and eigenvectors of the
+matrix A. The out parameter Values is the same as that obtained by calling the
+function Values. The out parameter Vectors is a matrix whose columns are the
+eigenvectors of the matrix A. The order of the columns corresponds to the
+order of the eigenvalues. The exception Constraint_Error is
+raised if A'Length(1) is not equal to A'Length(2). The index ranges of the
+parameter Vectors are those of A.
+
+Implementation Requirements
+
+The accuracy of these subprograms is implementation defined.
+
+
+Implementation Permissions
+
+The nongeneric equivalent packages may, but need not, be actual
+instantiations of the generic package for the appropriate predefined type.
 
-To be concluded
 
 !example
 
-To be done
+To be done by someone else.
 
 !discussion
 
-To be done
+The packages closely follow those in 13813. However, the discussion has been
+considerably simplified by assuming the properties of the corresponding scalar
+operations such as those in Numerics.Complex_Types. This removes a lot of
+explicit mention of raising exceptions and accuracy. Also remarks that these
+are standard mathematical operations have been deleted when the meaning is
+given by other words in the description.
+
+The component by component operations of * / and ** on vectors have been
+deleted on the grounds that they are not useful. (They might be useful for
+manipulating arrays in general but we are concerned with arrays used as vectors
+for linear algebra.)
+
+Operations for vector products have been added, However, this might be
+considered curious because they only apply in three-dimensional space.
+
+?? And I put the inner outer and vector products together which might not have
+been a good idea.
+
+It is hoped that there is not too much confusion between component when
+applied to the parts of a complex number and component when applied to a part
+of an array. 13813 uses element for the latter but the proper Ada term for an
+element of an array or field of a record is of course component.
+
+Observe that the index range of the various arrays is Integer (rather than
+Natural or Positive). This permits negative indices and in particular arrays
+with symmetric index ranges about zero.
+
+The function Identity_Matrix of 13813 has been changed to Unit_Matrix on the
+grounds that the prime concern is with linear algebra and not group theory.
+
+Functions have been added to Numerics.Generic_Real_Arrays for matrix inversion
+and related operations. On the grounds that it would be helpful to have simple
+algorithms available to users rather than none at all, no accuracy
+requirements are specified. Instead the accuracy is stated to be
+implementation-defined which means that the implementation must document it.
+
+The names chosen are Solve and Inverse. The former is overloaded, one version
+solves a single set of linear equations; the second solves multiple sets. Note
+that Inverse is not strictly necesssary becuase the effect can be obtained by
+using Solve and a Unit_Matrix thus
+
+I := Unit_Matrix(A'Length); B := Solve(A, I); -- same as B := Inverse (A);
+
+A function for computing the determinant has been added since it is helpful in
+deciding whether an array is ill-conditioned and therefore the results of
+inversion might be suspect.
+
+Similar functions have also been added for complex arrays.
+
+In addition, subprograms have been added for the computation of eigenvalues
+and vectors of real and complex arrays. Because these are somewhat more
+specialized than the elementary operations they have been placed in generic
+child packages such as Numerics.Generic_Real_Arrays.Eigen. The individual
+sunprograms are Values and Values_And_Vectors which means that they can be
+invoked by statements such as
+
+V := Eigen.Values(A);
+
+There is no separate function for eigenvectors since it is unlikely that these
+would be required without the eigenvalues.
+
+The most common application is with real symmetric matrices and these have real
+eigenvalues and vectors. Hence the subprograms in the real package only cater
+for symmetric matrices. Applications include: moments of inertia, elasticity,
+confidence regions and so on.
+
+A general real matrix can have complex eigenvalues and therefore provision for
+these and complex matrices is in the complex package. Note that the analogue of
+symmetric matrices in the complex domain is hermitian matrices (the transpose
+of a hermitian matrix equals its complex conjugate); these have real
+eigenvalues and distinct provision is made for these. The various overloadings
+are conveniently distinguished by the parameter profiles. (but could not be if
+there were separate fuinctions for the eigenvectors).
+
+?? I am somewhat rusty on eigensystems of non-symmetric matrices. I have said
+nothing about the normalization of the vectors. Note that I put all the eigen
+stuff in 3.3 - that's partly so it can be removed easily!!
+
+The original text of 13813 contains the remark "no complex conjugation is
+performed" in the description of multiplication of complex vectors to produce
+a complex matrix. I removed it because it seemed unhelpful to the user. The
+origin of this according to Donald Sando is that the NUMWG considered including
+various combined operations involving both multiplication and conjugation and
+indeed multiplication and transposition. Multiplication and conjugation are
+natural when dealing with for example Hermitian matrices. however, to provided
+all the appropriate combinations would have greatly increased the size of the
+package so it says in the rationale section F.3 of the 1998 version of 13813.
+
+Terry Froggatt adds further thoughts on this matter (private communication,
+17 July 2001. Thus
+
+2. THOUGHTS ON TRANSPOSITION.
+
+This comment arises from the last paragraph of F.3 of the Rationale in
+Annexe F ....
+
+In the Harrier flight-program we have several spatial co-ordinate systems,
+for example:
+
+Earth Axes: North - East - Up.
+Heading - Cross-heading - Up, (yawed wrt Earth).
+Aircraft Axes (rolled/pitched/yawed wrt Earth).
+HUD Axes (6 degrees down from Aircraft Axes).
+
+Conversion of vectors between pairs of these systems uses matrix
+multiplication one way, and multiplication by the tranpose of the same
+matrix the other way.
+
+Whilst X*CONJUGATE(Y) may be acceptable, since X and Y are the same size,
+the costs for X*TRANSPOSE(Y) are different: Y is generally much bigger than
+X, and it seems rather wasteful (mainly in time) to put the transpose of a
+matrix onto the stack then throw it away.
+
+An alternative is to calculate T := TRANSPOSE(Y) as soon as Y is calculated,
+so that it is available to convert several vectors between a pair of axis
+systems. This is wasteful (mainly in space) because several transposes have
+to be kept around simultaneously.
+
+What we actually use is one routine which implements X*TRANSPOSE(Y) directly:
+it performs the matrix multiplication but traverses Y in transpose-order. This
+takes no more time than the normal multiply, and no space for transposes (just
+code space).
+
+This is exactly what 13813's rationale rejected.
+
+Note that it should be possible to express the Harrier routine, and calls to
+it, quite neatly in Ada. Following the style of "Imaginary", define:
+
+type Transposed_Real_Matrix is private;
+...
+type Transposed_Real_Matrix is new Real_Matrix;
+function Transpose_In_Situ is new Unchecked_Conversion
+(Real_Matrix, Transposed_Real_Matrix);
+
+function "*" (Left: Real_Vector; Right: Transposed_Real_Matrix)
+                                        return (Real_Vector);
+
+This "*" knows it has to traverse Right in transposed order, and can be
+called (using operator notation, and no qualification) by
+
+"Z := X*Transpose_In_Situ(Y);".
+
+END of Remark by Terry Froggatt.
+
+Well now. It seems to me that playing with Hermitian matrices is pretty
+rare and so we can forget about possible combinations of conjugate and
+multiply. However, maybe there is a case for adding some combined
+multiplication and transposition for the real matrices as explained above.
 
 
 !ACATS Test

Questions? Ask the ACAA Technical Agent