CVS difference for ais/ai-00296.txt

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

--- ais/ai-00296.txt	2003/07/24 22:51:57	1.7
+++ ais/ai-00296.txt	2003/07/24 23:05:47	1.8
@@ -1,4 +1,4 @@
-!standard G (00)                                        03-07-22  AI95-00296/03
+!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
@@ -35,8 +35,9 @@
 renumbering of other clauses.
 
 In addition to the main facilities of 13813, these packages also include
-operations for matrix inversion and determinants. Furthermore, child packages
-provide subprograms for the determination of eigenvalues and eigenvectors.
+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
 
@@ -70,6 +71,7 @@
 
    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
 
@@ -120,12 +122,24 @@
    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
 
@@ -133,6 +147,10 @@
                          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
@@ -142,12 +160,15 @@
 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
+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.
+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
@@ -269,33 +290,65 @@
 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;
 
-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
+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;
 
-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
+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;
 
-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.
+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;
 
-This function returns the determinant of the matrix A. Constraint_Error is
-raised if A'Length(1) is not equal to A'Length(2).
+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;
@@ -308,8 +361,8 @@
 
 Implementation Requirements
 
-Accuracy requirements for the functions Solve, Inverse and Determinant are
-implementation defined.
+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
@@ -330,8 +383,14 @@
 
 Implementation Advice
 
-Implementations should implement the Solve and Inverse functions using
-established techniques such as Gaussian elimination with row interchanges.
+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
@@ -360,6 +419,7 @@
    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
 
@@ -528,15 +588,25 @@
    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;
 
-??  I have not added mixed real and complex for Solve.
-
    -- Other Complex_Matrix operations
 
    function Unit_Matrix (Order            : Positive;
@@ -550,19 +620,21 @@
 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.
+Complex type exported by Numerics.Complex_Types is systematically
+substituted for Complex, throughout. Nongeneric equivalents for each of
+the other predefined floating point types are defined similarly, with the
+names Numerics.Short_Complex_Arrays, Numerics.Long_Complex_Arrays, etc.
 
-Two types are defined and exported by Ada.Numerics.Generic_Complex_Arrays.
+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.
+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
@@ -907,33 +979,65 @@
 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;
 
-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
+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;
 
-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
+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;
 
-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.
+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;
 
-This function returns the determinant of the matrix A. Constraint_Error is
-raised if A'Length(1) is not equal to A'Length(2).
+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)
@@ -947,8 +1051,8 @@
 
 Implementation Requirements
 
-Accuracy requirements for the functions Solve, Inverse and Determinant are
-implementation defined.
+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
@@ -973,8 +1077,14 @@
 
 Implementation Advice
 
-Implementations should implement the Solve and Inverse functions using
-established techniques such as Gaussian elimination with row interchanges.
+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
@@ -997,10 +1107,14 @@
    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;
 
@@ -1014,27 +1128,32 @@
 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;
 
-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).
+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);
 
-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
+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.
-
-For both subprograms, if the matrix A is not symmetric then Argument_Error is
-raised.
+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
@@ -1047,26 +1166,38 @@
    -- 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;
 
@@ -1081,70 +1212,88 @@
 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;
 
-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).
+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);
 
-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
+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;
 
-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 Argument_Error is raised if the matrix
-is not Hermitian.
+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);
 
-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
+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;
 
-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).
+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);
 
-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.
+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
@@ -1162,10 +1311,44 @@
 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.
 
-To be done by someone else.
+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
 
@@ -1213,6 +1396,26 @@
 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.
 
@@ -1245,18 +1448,16 @@
 Hermitian no requirement has been placed on the independence or modulus of
 the vectors.
 
-Consideration was given to the inclusion of explicit facilities for LU
-decomposition (as provided for example in the so-called Basic Linear Algebra
-Subprograms (BLAS)). This permits more rapid computation when several
-operations such as solving equations, determinant evaluation and eigensystems
-are to be performed using the same matrix. However, modern hardware is so
-fast that this would only seem to be necessary for very large sets of equations
-which are not the target of these simple facilities. Moroever, adding LU
-decomposition introduces complexity for the user.
-
-Consideration was also given to a fuller implementation of the BLAS. However,
-this seems out of place in a language standard since it would be extremely long
-and specialized. Such a fuller interface to the BLAS could be provided in
+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
@@ -1604,6 +1805,13 @@
 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.
 
 ****************************************************************
 

Questions? Ask the ACAA Technical Agent