--- ais/ai-00296.txt 2003/07/26 03:26:01 1.9 +++ ais/ai-00296.txt 2003/09/19 01:42:27 1.10 @@ -1,4 +1,4 @@ -!standard G (00) 03-07-23 AI95-00296/04 +!standard G (00) 03-09-01 AI95-00296/05 !class amendment 02-06-07 !status work item 03-01-23 !status received 02-06-07 @@ -35,9 +35,10 @@ renumbering of other clauses. In addition to the main facilities of 13813, these packages also include -operations for matrix inversion and determinants including the use of -LU decomposition. Furthermore, child packages provide subprograms for the -determination of eigenvalues and eigenvectors. +operations for the solution of linear equations, matrix inversion and +determinants. Furthermore, child packages provide subprograms for the +determination of the eigenvalues and eigenvectors of real symmetric matrices +and Hermitian matrices. !wording @@ -53,8 +54,9 @@ equivalents of these packages for each of the predefined floating point types are also provided as children of Numerics. -Generic child packages for the determination of eigenvalues and eigenvectors -are defined in G.3.3. +Generic child packages for the determination of the eigenvalues and +eigenvectors of real symmetric matrices and Hermitian matrices are defined +in G.3.3. G.3.1 Real Vectors and Matrices @@ -71,7 +73,6 @@ 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 @@ -122,24 +123,12 @@ 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 @@ -160,15 +149,12 @@ point types are defined similarly, with the names Numerics.Short_Real_Arrays, Numerics.Long_Real_Arrays, etc. -Three types are defined and exported by Ada.Numerics.Generic_Real_Arrays. The +Two types are defined and exported by Ada.Numerics.Generic_Real_Arrays. The composite type Real_Vector is provided to represent a vector with components of type Real; it is defined as an unconstrained, one-dimensional array with an index of type Integer. The composite type Real_Matrix is provided to represent a matrix with components of type Real; it is defined as an unconstrained, -two-dimensional array with indices of type Integer. The 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. +two-dimensional array with indices of type Integer. The effect of the various functions is as described below. In most cases the functions are described in terms of corresponding scalar operations of the type @@ -290,65 +276,33 @@ the type Real to each component of Left and to the scalar Right. The index ranges of the matrix result are those of Left. -function Decompose (A : Real_Matrix) return Real_LU; - -This function returns the LU decomposed form of the matrix A. Constraint_Error -is raised if A'Length(1) and A'Length(2) are not equal. - -function Lower (LU: Real_LU) return Real_Matrix; - -This function returns the lower triangular matrix of the decomposed matrix LU. -The diagonal components are 1.0. Both index ranges of the matrix result are -A'Range(1) where A is the matrix of which LU is the decomposition. - -function Upper (LU: Real_LU) return Real_Matrix; - -This function returns the upper triangular matrix of the decomposed matrix LU. -The index ranges of the matrix result are those of A where A is the matrix of -which LU is the decomposition. - -function Permutation (LU: Real_LU) return Real_Matrix; - -This function returns a matrix whose components are zero except for one -component in each row and column which is 1.0. It represents the row -interchanges performed during the decomposition. Both index ranges of the -result are A'Range(1) where A is the matrix of which LU is the decomposition. - function Solve (A : Real_Matrix; X: Real_Vector) return Real_Vector; -function Solve (LU : Real_LU; X: Real_Vector) return Real_Vector; -These functions return a vector Y such that X is (nearly) equal to A * Y where A -is either provided directly or is the matrix of which LU is the decomposition. -This is the standard mathematical operation for solving a single set of linear +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; -function Solve (LU : Real_LU; X : Real_Matrix) return Real_Matrix; -These functions return a matrix Y such that X is (nearly) equal to A * Y where A -is either provided directly or is the matrix of which LU is the decomposition. -This is the standard mathematical operation for solving several sets of linear +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; -function Inverse (LU : Real_LU) return Real_Matrix; -These functions return a matrix B such that A * B is (nearly) the unit matrix -where A is either provided directly or is the matrix of which LU is the -decomposition. The index ranges of the result are those of A. Constraint_Error -is raised if A'Length(1) is not equal to A'Length(2). Constraint_Error is -raised if the matrix A is ill-conditioned. +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; -function Determinant (LU : Real_LU) return Real'Base; -These functions return the determinant of the matrix A where A is either -provided directly or is the matrix of which LU is the decomposition. -Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). +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 Real_Matrix; @@ -361,8 +315,8 @@ Implementation Requirements -Accuracy requirements for the functions Decompose, Lower, Upper, Solve, -Inverse and Determinant are implementation defined. +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 @@ -383,14 +337,11 @@ Implementation Advice -Implementations should implement the Solve and Inverse functions using LU -decomposition with row interchanges. This could be via a binding to an -implementation of the so-called BLAS (Basic Linear Algebra Subprograms). -Bindings to other facilities in the BLAS could be provided in child packages. - -The implementation of the type Real_LU must retain information regarding the -ranges of the original array since these are needed by Solve and Inverse for -compatibility with the direct versions. +Implementations should implement the Solve and Inverse functions using +established techniques such as LU decomposition with row interchanges followed +by back and forward substitution. Implementations may choose to refine the +result by performing an iteration on the residuals; if this is done then it +should be documented. It is not the intention that any special provision should be made to determine whether a matrix is ill-conditioned or not. The naturally occurring @@ -419,7 +370,6 @@ 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 @@ -588,24 +538,12 @@ function "/" (Left : Complex_Matrix; Right : Real'Base) return Complex_Matrix; - -- Complex_Matrix LU decomposition and related subprograms - - function Decompose (A : Complex_Matrix) return Complex_LU; - function Lower (LU: Complex_LU) return Complex_Matrix; - function Upper (LU: Complex_LU) return Complex_Matrix; - function Permutation (LU: Complex_LU) return Complex_Matrix; - -- Complex_Matrix inversion and related operations function Solve (A : Complex_Matrix; X: Complex_Vector) return Complex_Vector; - function Solve (LU : Complex_LU; X: Complex_Vector) return Complex_Vector; function Solve (A, X : Complex_Matrix) return Complex_Matrix; - function Solve (LU : Complex_LU; X : Complex_Matrix) return Complex_Matrix; - function Inverse (A : Complex_Matrix) return Complex_Matrix; - function Inverse (A : Complex_LU) return Complex_Matrix; function Determinant (A : Complex_Matrix) return Complex; - function Determinant (A : Complex_LU) return Complex; -- Other Complex_Matrix operations @@ -625,16 +563,13 @@ the other predefined floating point types are defined similarly, with the names Numerics.Short_Complex_Arrays, Numerics.Long_Complex_Arrays, etc. -Three types are defined and exported by Ada.Numerics.Generic_Complex_Arrays. +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 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. +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 @@ -979,65 +914,33 @@ Numerics.Generic_Complex_Types to each component of the matrix Left and the scalar Right. The index ranges of the result are those of Left. -function Decompose (A : Complex_Matrix) return Complex_LU; - -This function returns the LU decomposed form of the matrix A. Constraint_Error -is raised if A'Length(1) and A'Length(2) are not equal. - -function Lower (LU: Complex_LU) return Complex_Matrix; - -This function returns the lower triangular matrix of the decomposed matrix LU. -The diagonal components are 1.0. Both index ranges of the matrix result are -A'Range(1) where A is the matrix of which LU is the decomposition. - -function Upper (LU: Complex_LU) return Complex_Matrix; - -This function returns the upper triangular matrix of the decomposed matrix LU. -The index ranges of the matrix result are those of A where A is the matrix of -which LU is the decomposition. - -function Permutation (LU: Complex_LU) return Real_Matrix; - -This function returns a real matrix whose components are zero except for one -component in each row and column which is 1.0. It represents the row -interchanges performed during the decomposition. Both index ranges of the -result are A'Range(1) where A is the matrix of which LU is the decomposition. - function Solve (A : Complex_Matrix; X: Complex_Vector) return Complex_Vector; -function Solve (LU : Complex_LU; X: Complex_Vector) return Complex_Vector; -These functions return a vector Y such that X is (nearly) equal to A * Y where A -is either provided directly or is the matrix of which LU is the decomposition. -This is the standard mathematical operation for solving a single set of linear +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; -function Solve (LU : Complex_LU; X : Complex_Matrix) return Complex_Matrix; -These functions return a matrix Y such that X is (nearly) equal to A * Y where A -is either provided directly or is the matrix of which LU is the decomposition. -This is the standard mathematical operation for solving several sets of linear +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; -function Inverse (LU : Complex_LU) return Complex_Matrix; -These functions return a matrix B such that A * B is (nearly) the unit matrix -where A is either provided directly or is the matrix of which LU is the -decomposition. The index ranges of the result are those of A. Constraint_Error -is raised if A'Length(1) is not equal to A'Length(2). Constraint_Error is raised -if the matrix A is ill-conditioned. +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; -function Determinant (LU : Complex_LU) return Complex; -These functions return the determinant of the matrix A where A is either -provided directly or is the matrix of which LU is the decomposition. -Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). +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) @@ -1051,8 +954,8 @@ Implementation Requirements -Accuracy requirements for the functions Decompose, Lower, Upper, Solve, -Inverse and Determinant are implementation defined. +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 @@ -1077,14 +980,10 @@ Implementation Advice -Implementations should implement the Solve and Inverse functions using LU -decomposition with row interchanges. This could be via a binding to an -implementation of the so-called BLAS (Basic Linear Algebra Subprograms). -Bindings to other facilities in the BLAS could be provided in child packages. - -The implementation of the type Complex_LU must retain information regarding the -ranges of the original array since these are needed by Solve and Inverse for -compatibility with the direct versions. +Implementations should implement the Solve and Inverse functions using +established techniques. Implementations may choose to refine the result by +performing an iteration on the residuals; if this is done then it should be +documented. It is not the intention that any special provision should be made to determine whether a matrix is ill-conditioned or not. The naturally occurring @@ -1106,15 +1005,13 @@ package Ada.Numerics.Generic_Real_Arrays.Eigen is pragma Pure(Eigen); + -- Eigensystem of a real symmetric matrix + 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; @@ -1128,77 +1025,44 @@ The effect of the subprograms is as described below. function Values(A : Real_Matrix) return Real_Vector; -function Values(LU : Real_LU) return Real_Vector; -These functions return the eigenvalues of the symmetric matrix A as a vector -sorted into order with the largest first. The matrix A is either provided -directly or is the matrix of which LU is the decomposition. The exception -Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). The -index range of the result is A'Range(1). If the matrix A is not symmetric -then Argument_Error is raised. +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); -procedure Values_And_Vectors(LU : in Real_LU; - Values : out Real_Vector; - Vectors : out Real_Matrix); -These procedures compute both the eigenvalues and eigenvectors of the symmetric -matrix A. The matrix A is either provided directly or is the matrix of which LU -is the decomposition. The out parameter Values is the same as that obtained by -calling the function Values. The out parameter Vectors is a matrix whose columns -are the eigenvectors of the matrix A. The order of the columns corresponds to -the order of the eigenvalues. The eigenvectors are normalized and mutually +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. -If the matrix A is not symmetric then Argument_Error is raised. + +For both subprograms, if the matrix A is not symmetric then Argument_Error is +raised. -Eigenvalues and eigenvectors of a general real or complex matrix are computed -using subprograms in the child package -Ada.Numerics.Generic_Complex_Arrays.Eigen whose declaration is: +Eigenvalues and eigenvectors of a Hermitian matrix are computed using +subprograms in the child package Ada.Numerics.Generic_Complex_Arrays.Eigen +whose declaration is: generic package Ada.Numerics.Generic_Complex_Arrays.Eigen is pragma Pure(Eigen); - -- Eigensystem of a general real matrix - - function Values(A : Real_Matrix) return Complex_Vector; - function Values(LU : Real_LU) return Complex_Vector; - - procedure Values_And_Vectors(A : in Real_Matrix; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); - procedure Values_And_Vectors(LU : in Real_LU; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); - -- Eigensystem of a Hermitian matrix function Values(A : Complex_Matrix) return Real_Vector; - function Values(LU : Complex_LU) return Real_Vector; procedure Values_And_Vectors(A : in Complex_Matrix; Values : out Real_Vector; Vectors : out Complex_Matrix); - procedure Values_And_Vectors(LU : in Complex_LU; - Values : out Real_Vector; - Vectors : out Complex_Matrix); - -- Eigensystem of a general complex matrix - - function Values(A : Complex_Matrix) return Complex_Vector; - function Values(LU : Complex_LU) return Complex_Vector; - - procedure Values_And_Vectors(A : in Complex_Matrix; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); - procedure Values_And_Vectors(LU : in Complex_LU; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); - end Ada.Numerics.Generic_Complex_Arrays.Eigen; The library package Numerics.Complex_Arrays.Eigen is declared pure and defines @@ -1211,84 +1075,28 @@ The effect of the subprograms is as described below. -function Values(A : Real_Matrix) return Complex_Vector; -function Values(LU : Real_LU) return Complex_Vector; - -These functions return the eigenvalues of the matrix A as a vector sorted by -order of modulus with the largest first. The matrix A is either provided -directly or is the matrix of which LU is the decomposition. The exception -Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). The -index range of the result is A'Range(1). - -procedure Values_And_Vectors(A : in Real_Matrix; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); -procedure Values_And_Vectors(LU : in Real_LU; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); - -These procedures compute both the eigenvalues and eigenvectors of the -matrix A. The matrix A is either provided directly or is the matrix of which LU -is the decomposition. The out parameter Values is the same as that obtained by -calling the function Values. The out parameter Vectors is a matrix whose columns -are the eigenvectors of the matrix A. The order of the columns corresponds to -the order of the eigenvalues. The exception Constraint_Error is -raised if A'Length(1) is not equal to A'Length(2). The index ranges of the -parameter Vectors are those of A. - function Values(A : Complex_Matrix) return Real_Vector; -function Values(LU : Complex_LU) return Real_Vector; -These functions return the eigenvalues of the Hermitian matrix A as a vector -sorted into order with the largest first. The matrix A is either provided -directly or is the matrix of which LU is the decomposition. The exception -Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). The -index range of the result is A'Range(1). The exception Argument_Error is -raised if the matrix is not Hermitian. +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). procedure Values_And_Vectors(A : in Complex_Matrix; Values : out Real_Vector; Vectors : out Complex_Matrix); -procedure Values_And_Vectors(LU : in Complex_LU; - Values : out Real_Vector; - Vectors : out Complex_Matrix); -These procedures compute both the eigenvalues and eigenvectors of the Hermitian -matrix A. The matrix A is either provided directly or is the matrix of which LU -is the decomposition. The out parameter Values is the same as that obtained by -calling the function Values. The out parameter Vectors is a matrix whose columns -are the eigenvectors of the matrix A. The order of the columns corresponds to -the order of the eigenvalues. The eigenvectors are mutually orthonormal, +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 Argument_Error is raised if -the matrix is not Hermitian. - -function Values(A : Complex_Matrix) return Complex_Vector; -function Values(LU : Complex_LU) return Complex_Vector; - -These functions return the eigenvalues of the matrix A as a vector sorted by -order of modulus with the largest first. The matrix A is either provided -directly or is the matrix of which LU is the decomposition. The exception -Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). The index -range of the result is A'Range(1). - -procedure Values_And_Vectors(A : in Complex_Matrix; - Values : out Complex_Vector; - - Vectors : out Complex_Matrix); -procedure Values_And_Vectors(LU : in Complex_LU; - Values : out Complex_Vector; - Vectors : out Complex_Matrix); +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. +For both subprograms, if the matrix A is not Hermitian then Argument_Error is +raised. Implementation Requirements @@ -1311,47 +1119,32 @@ negation is exact. If it doesn't then the paragraph above needs changing or G.2.1 needs changing. -AARM Note - -In this child package, Generic_Complex_Arrays.Eigen, the full view of the type -Real_LU is not available. However, the details may be obtained by calling the -functions Lower, Upper and Permutation in the formal package Real_Arrays. This -is the real reason for declaring these subprograms although they might be -helpful to the inquisitive user. - !example -This example illustrates the use of the LU decomposed form to avoid -unnecessary computation. -with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; -with Ada.Numerics.Real_Arrays.Eigen; -procedure An_Example is - N : Integer := ... - Mat : Real_Array(1 .. N, 1 .. N); - P, Q, R : Real_Vector(1 .. N); - X, Y, Z : Real_Vector(1 .. N); - -- compute some value for Mat and values for vectors P, Q and R; - - begin - LUMat : constant Real_LU := Decompose(Mat); - -- Now compute various properties of Mat - D : Float; - Inverse_Mat : Real_Array(1 .. N); - Mat_Eigenvalues : Real_Array(1 .. N); - ... - D := Determinant(LUMat); - Inverse_Mat := Invert(LUMat); - X := Solve(LUMat, P); - Y := Solve(LUMat, Q); - Z := Solve(LUMat, R); - Mat_Eigenvalues := Eigen.Values(LUMat); - ... - end; -end An_Example; - !discussion +Section G.1.1 of the Rationale for Ada 95 says + +"A decision was made to abbreviate the Ada 95 packages by omitting the vector +and matrix types and operations. One reason was that such types and operations +were largely self-evident, so that little real help would be provided by +defining them in the language. Another reason was that a future version of Ada +might add enhancements for array manipulation and so it would be inappropriate +to lock in such operations prematurely." + +It is now clear that such enhancements will not be added so the second +justification for omitting the facilities of 13813 disappears. In order to +overcome the objection that everything is self-evident we have taken the +approach that we should add some basic facilities that are commonly required, +not completely trivial to implement but on the other hand are mathematically +well understood. + +The overall goal is thus to provide something useful for the user who is not a +numerical professional and also to provide a baseline of types and operations +that form a firm foundation for binding to more general facilities such as the +well-known BLAS (Basic Linear Algebra Subprograms, see www.netlib.org/blas). + The packages closely follow those in 13813. However, the discussion has been considerably simplified by assuming the properties of the corresponding scalar operations such as those in Numerics.Complex_Types. This removes a lot of @@ -1379,6 +1172,15 @@ The function Identity_Matrix of 13813 has been changed to Unit_Matrix on the grounds that the prime concern is with linear algebra and not group theory. +The accuracy of most simple operations follows from the underlying operations +on scalar components. In the case of inner product there is the potential for +special operations to improve the speed and/or accuracy. We have simply +required that in the exact mode, the accuracy should be no less than that of +the canonical implementation using a loop statement. This is because on some +hardare, built-in instructions which are fast actually lose accuracy. Note +that the Fortran language standard recognizes the existence of inner product +as a special case but imposes no accuracy requirements at all. + 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 @@ -1392,39 +1194,41 @@ I := Unit_Matrix(A'Length); B := Solve(A, I); -- same as B := Inverse (A); +A common technique when solving sets of linear equations is to refine the result +by an iteration on the residuals. Thus to solve the set of equations Ax = y we +first perform + +X := Solve(A, Y); + +we then compute the error D by multiplying back thus + +D := Y - A*X; + +and then compute the refinement to X + +DX := Solve(A, D); + +and then correct X + +X := X + DX; + +Implementations may choose to do this automatically; it requires little extra +computation if LU decomposition is used internally. If they do such refinement +then it should be documented. + A function for computing the determinant has been added since it is helpful in deciding whether an array is ill-conditioned and therefore the results of inversion might be suspect. -In addition a private type has been added to hold the LU decomposed form of a -matrix. This enables repeated calls of Solve etc to be done efficiently by -avoiding unnecessary repeated decomposition. Further overloadings of the -various subprograms using the LU form are thus provided. This private type has -unknown discriminants and thus objects of this type must be initialized with a -call of Decompose. The use of a private type prevents the user from making -errors by confusing the original and LU forms of a matrix. - -The private type has unknown discriminants. The full type must retain the -information regarding the bounds of the original array since it is needed for -some subprograms. The functions Lower, Upper and Permutation give access to the -details of the decomposition. Note that the permutation describing the row -interchanges is given in the rather unorthodox form of a real matrix. It would -be more conventional just to provide an integer vector defining the permutation -but there is no convenient type available and declaring an integer vector type -just for this seems tiresome. Moreover, the matrix form of the permutation -enables the original array to be reconstructed simply by multiplication. - -?? I hope this is true. Someone better check this carefully. - Similar functions have also been added for complex arrays. However, it was not deemed necessary to provide for mixed real and complex operands for Solve. In addition, subprograms have been added for the computation of eigenvalues -and vectors of real and complex arrays. Because these are somewhat more -specialized than the elementary operations they have been placed in generic -child packages such as Numerics.Generic_Real_Arrays.Eigen. The individual -subprograms are Values and Values_And_Vectors which means that they can be -invoked by statements such as +and vectors of real symmetric matrices and Hermitian matrices. Because these +are somewhat more specialized than the elementary operations they have been +placed in generic child packages such as Numerics.Generic_Real_Arrays.Eigen. +The individual subprograms are Values and Values_And_Vectors which means that +they can be invoked by statements such as V := Eigen.Values(A); @@ -1432,36 +1236,83 @@ 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, +eigenvalues and vectors. The subprograms for these are in the real child +package. Applications include: moments of inertia, elasticity, confidence regions and so on. -A general real matrix can have complex eigenvalues and therefore provision for -these and complex matrices is in the complex package. Note that the analogue of -symmetric matrices in the complex domain is Hermitian matrices (the transpose -of a Hermitian matrix equals its complex conjugate); these have real -eigenvalues and distinct provision is made for these. The various overloadings -are conveniently distinguished by the parameter profiles (but could not be if -there were separate functions for the eigenvectors). - -In the case of eigenvectors of matrices which are neither symmetric nor -Hermitian no requirement has been placed on the independence or modulus of -the vectors. - -All the eigenvalue and vector subprograms also have overloadings for using the -LU decomposed forms. A curiosity arises in the case of the eigenvalues of a -general real matrix because this needs a parameter of type Real_LU. However, -the full view of this type is not available in the complex eigenvalue package -and so the decomposed form has to be accessed via the subprograms Lower, Upper -and Permutation. - -Consideration was given to a fuller implementation of the BLAS. However, this -seems out of place in a language standard since it would be extremely long and -specialized. Such a fuller interface to the BLAS could be provided in +A slight problem arises in deciding who should check that a matrix is symmetric, +the user or the package. Computations of, for example, a covariance matrix will +result in a matrix that ought to be exactly symmetric but small errors might be +introduced which mean that it is not exactly symmetric. The onus could be +placed on the user to ensure that the matrix is exactly symmetric or +alternatively some tolerance could be passed to the eigen subprograms so that +they can perform a reasonable check. Passing a tolerance level adds an +irritating parameter, raises the issue of how the tolerance should be expressed +and gives the user the problem of what tolerance to request. Moreover, the +algorithms still have to decide which actual values should be used for those +pairs of components that are not exactly equal. + +Accordingly, we have placed the onus on the user to ensure that the matrix is +exactly symmetric. This could be done for example by taking the mean of the +matrix and its transpose. The test in the subprograms can then test for exact +equality which is guaranteed to deliver the correct answer. If the test fails +then Argument_Error is raised. This exception, declared in Ada.Numerics, is +the normal exception for numeric operations when the argument is out of range. +Note that errors such as when bounds of arrays do not match raise +Constraint_Error by analogy with built-in array operations. + +The complex equivalent of real symmetric matrices are Hermitian matrices. +Hermitian matrices are such that their transpose (that is with rows and columns +interchanged) equals their complex conjugate (that is with the sign of imaginary +parts reversed). Hermitian matrices also have real eigenvalues and vectors. The +subprograms for these are in the complex child package. Applications include +Quantum Mechanics. + +Again we have placed the onus on the user to ensure that the matrix is +Hermitian. The check in the package can then use strict equality for the real +parts and negation followed by equality for the imaginary parts. + +We considered providing subprograms for the determination of eigenvalues and +eigenvectors of general real and complex matrices. Such matrices can have +complex eigenvalues and therefore provision for these would have to be in the +complex package. However, there are mathematical difficulties with these general +cases which are in strong contrast to the real symmetric and Hermitian matrices. +Thus, Numerical Recipes by Press, Flannery, Teukolsky and Vetterling says +regarding the real case: + +"The algorithms for symmetric matrices ... are highly satisfactory in practice. +By contrast, it is impossible to design equally satisfactory algorithms for the +nonsymmetric case. There are two reasons for this. First, the eigenvalues of a +nonsymmetric matrix can be very sensitive to small changes in the matrix +elements. Second, the matrix itself can be defective so that there is no +complete set of eigenvectors. We emphasize that these difficulties are intrinsic +properties of certain nonsymmetric matrices, and no numerical procedure can cure +them." + +Similar remarks apply to complex matrices where Hermitian matrices are well- +behaved but non-Hermitian matrices can be troublesome. + +In view of these computational difficulties and the fact that requiring the +eigensystem of general matrices is uncommon, we decided not to provide such +facilities. + +Consideration was also given to the inclusion of explicit facilities for LU +decomposition (as provided for example in the BLAS). LU decomposition is a +common first step for many operations. Thus making it available to the user +permits more rapid computation when several operations such as solving +equations, 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 and these are not the +target of these simple facilities. Moreover, adding explicit LU decomposition +introduces complexity for the user. + +Consideration was also given to a fuller implementation of the BLAS. However, +this seems out of place in a language standard since it would be extremely long +and specialized. Such a fuller interface to the BLAS could be provided in additional child packages. The goal here has been to provide convenient access to simple type declarations and the more commonly required operations that are -not easy for the user to program. These operations can of course be implemented -as a binding to an implementation of part of the BLAS. +not trivial for the user to program. These operations can of course be +implemented as a binding to an implementation of part of the BLAS. !ACATS Test @@ -1889,6 +1740,51 @@ We can discuss it in Sydney. But any thoughts before then might be welcome. I suppose I will have to do the wording for downward closures now - groan. + +**************************************************************** + +From: John Barnes +Sent: Tuesday, September 2, 2003 6:40 AM + +This version [version /05 - ED] 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. Moreover, the eigenvalues and vectors of +nonsymmetric, non-Hermitian matrices have been removed because of potential +computational difficulties. + +**************************************************************** + +From: John Barnes +Sent: Tuesday, September 2, 2003 6:54 AM + +I have had a number of further discussions on this topic +with people at NAG, Farnborough, Pascal etc. I have also +invested in a heavy tome on numerical stuff. + +As a consequence I have concluded that making LU +decomposition visible to the user was unnecessary. Moreover, +after reading around, it seems that some of the eigenvalue +calculations can be tricky and accordingly I have reduced +the complex package to the one straightforward case. + +I attach a paper written for the Ada User Journal which +explains this in more detail. + +And so I have written yet another version of this AI with +which I feel fairly comfortable. This is version 5 dated +2003/09/01 and I have sent this to Randy for placing on the +database. For those who cannot wait I also attach it. + +For those of a historical inclination, the versions are as +follows + +02 as presented at Toulouse + +03 as 02 plus some minor corrections agreed at Toulouse + +04 as 03 plus LU decomposition + +05 as 03 minus eigensystems of general matrices ****************************************************************

Questions? Ask the ACAA Technical Agent