--- ais/ai-00296.txt 2003/09/25 01:49:11 1.11 +++ ais/ai-00296.txt 2003/10/29 22:54:11 1.12 @@ -1,5 +1,9 @@ -!standard G (00) 03-09-01 AI95-00296/05 +!standard G.3 (01) 03-10-24 AI95-00296/06 +!standard G.3.1 (01) +!standard G.3.2 (01) !class amendment 02-06-07 +!status work item 03-10-29 +!status ARG Approved 10-0-0 03-10-03 !status work item 03-01-23 !status received 02-06-07 !priority Medium @@ -35,16 +39,15 @@ renumbering of other clauses. In addition to the main facilities of 13813, these packages also include -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. +subprograms for the solution of linear equations, matrix inversion, +determinants, and the determination of the eigenvalues and eigenvectors of +real symmetric matrices and Hermitian matrices. !wording Add a new clause G.3 as follows -G.3 Vector and matrix manipulation +G.3 Vector and Matrix Manipulation Types and operations for the manipulation of real vectors and matrices are provided in Generic_Real_Arrays, which is defined in G.3.1. Types and @@ -54,10 +57,6 @@ 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 the eigenvalues and -eigenvectors of real symmetric matrices and Hermitian matrices are defined -in G.3.3. - G.3.1 Real Vectors and Matrices The generic library package Numerics.Generic_Real_Arrays has the following @@ -130,6 +129,14 @@ function Inverse (A : Real_Matrix) return Real_Matrix; function Determinant (A : Real_Matrix) return Real'Base; + -- Eigenvalues and vectors of a real symmetric matrix + + function Eigenvalues(A : Real_Matrix) return Real_Vector; + + procedure Eigensystem(A : in Real_Matrix; + Values : out Real_Vector; + Vectors : out Real_Matrix); + -- Other Real_Matrix operations function Unit_Matrix (Order : Positive; @@ -304,6 +311,28 @@ This function returns the determinant of the matrix A. Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). +function Eigenvalues(A : Real_Matrix) return Real_Vector; + +This function returns the eigenvalues of the symmetric matrix A as a vector +sorted into order with the largest first. 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 A +is not symmetric. + +procedure Eigensystem(A : in Real_Matrix; + Values : out Real_Vector; + Vectors : out Real_Matrix); + +This procedure computes both the eigenvalues and eigenvectors of the symmetric +matrix A. The out parameter Values is the same as that obtained by calling the +function Eigenvalues. The out parameter Vectors is a matrix whose columns are +the eigenvectors of the matrix A. The order of the columns corresponds to the +order of the eigenvalues. The eigenvectors are normalized and mutually +orthogonal (they are orthonormal), including when there are repeated +eigenvalues. 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 A is not symmetric. + function Unit_Matrix (Order : Positive; First_1, First_2 : Integer := 1) return Real_Matrix; @@ -315,21 +344,33 @@ Implementation Requirements -Accuracy requirements for the functions Solve, Inverse and Determinant are -implementation defined. +Accuracy requirements for the subprograms Solve, Inverse, Determinant, +Eigenvalues and Eigensystem are implementation defined. For operations not involving an inner product, the accuracy requirements are those of the corresponding operations of the type Real in both the strict mode and the relaxed mode (see G.2). For operations involving an inner product, no requirements are specified in -the relaxed mode. In the strict mode the accuracy should be at least that of -the canonical implementation of multiplication and addition using the -corresponding operations of type Real'Base and performing the cumulative -addition using ascending indices. Implementations shall document any -techniques used to reduce cancellation errors such as extended precision -arithmetic. +the relaxed mode. In the strict mode the modulus of the absolute error of +the inner product X*Y shall not exceed g*abs(X)*abs(Y) where g is defined as + g = X'Length * Machine_Radix**(1-Machine_Mantissa) + +Implementations shall document any techniques used to reduce cancellation +errors such as extended precision arithmetic. + +AARM Note + +The above accuracy requirement is met by the canonical implementation of the +inner product by multiplication and addition using the corresponding +operations of type Real'Base and performing the cumulative addition using +ascending indices. Note however, that some hardware provides special +operations for the computation of the inner product and although these may be +fast they may not meet the accuracy requirement specified. See Accuracy and +Stability of Numerical Algorithms By N J Higham (ISBN 0-89871-355-2), +Section 3.1. + Implementation Permissions The nongeneric equivalent packages may, but need not, be actual @@ -339,7 +380,7 @@ 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 +by back and forward substitution. Implementations are recommended to refine the result by performing an iteration on the residuals; if this is done then it should be documented. @@ -349,6 +390,8 @@ functions with an ill-conditioned matrix and thus raise Constraint_Error is sufficient. +The test that a matrix is symmetric may be performed by using the equality +operator to compare the relevant components. G.3.2 Complex Vectors and Matrices @@ -545,6 +588,14 @@ function Inverse (A : Complex_Matrix) return Complex_Matrix; function Determinant (A : Complex_Matrix) return Complex; + -- Eigenvalues and vectors of a Hermitian matrix + + function Eigenvalues(A : Complex_Matrix) return Real_Vector; + + procedure Eigensystem(A : in Complex_Matrix; + Values : out Real_Vector; + Vectors : out Complex_Matrix); + -- Other Complex_Matrix operations function Unit_Matrix (Order : Positive; @@ -942,6 +993,28 @@ This function returns the determinant of the matrix A. Constraint_Error is raised if A'Length(1) is not equal to A'Length(2). +function Eigenvalues(A : Complex_Matrix) return Real_Vector; + +This function returns the eigenvalues of the Hermitian matrix A as a vector +sorted into order with the largest first. 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 A +is not Hermitian. + +procedure Eigensystem(A : in Complex_Matrix; + Values : out Real_Vector; + Vectors : out Complex_Matrix); + +This procedure computes both the eigenvalues and eigenvectors of the Hermitian +matrix A. The out parameter Values is the same as that obtained by calling the +function Eigenvalues. The out parameter Vectors is a matrix whose columns are +the eigenvectors of the matrix A. The order of the columns corresponds to the +order of the eigenvalues. The eigenvectors are mutually orthonormal, +including when there are repeated eigenvalues. 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 A is not Hermitian. + function Unit_Matrix (Order : Positive; First_1, First_2 : Integer := 1) return Complex_Matrix; @@ -954,21 +1027,37 @@ Implementation Requirements -Accuracy requirements for the functions Solve, Inverse and Determinant are -implementation defined. +Accuracy requirements for the subprograms Solve, Inverse, Determinant, +Eigenvalues and Eigensystem are implementation defined. For operations not involving an inner product, the accuracy requirements are those of the corresponding operations of the type Real'Base and Complex in both the strict mode and the relaxed mode (see G.2). For operations involving an inner product, no requirements are specified in -the relaxed mode. In the strict mode the accuracy should be at least that of -the canonical implementation of multiplication and addition using the -corresponding operations of type Real'Base and Complex and performing the -cumulative addition using ascending indices. Implementations shall document -any techniques used to reduce cancellation errors such as extended precision -arithmetic. +the relaxed mode. In the strict mode the modulus of the absolute error of the +inner product X*Y shall not exceed g*abs(X)*abs(Y) where g is defined as + + g = sqrt(2.0) * X'Length * Machine_Radix**(1-Machine_Mantissa) for mixed + complex and real operands + g = sqrt(8.0) * X'Length * Machine_Radix**(1-Machine_Mantissa) for two + complex operands + +Implementations shall document any techniques used to reduce cancellation +errors such as extended precision arithmetic. + +AARM Note + +The above accuracy requirement is met by the canonical implementation of the +inner product by multiplication and addition using the corresponding +operations of type Complex and performing the cumulative addition using +ascending indices. Note however, that some hardware provides special +operations for the computation of the inner product and although these may be +fast they may not meet the accuracy requirement specified. See Accuracy and +Stability of Numerical Algorithms by N J Higham (ISBN 0-89871-355-2), +Sections 3.1 and 3.6. + Implementation Permissions The nongeneric equivalent packages may, but need not, be actual @@ -981,7 +1070,7 @@ Implementation Advice Implementations should implement the Solve and Inverse functions using -established techniques. Implementations may choose to refine the result by +established techniques. Implementations are recommended to refine the result by performing an iteration on the residuals; if this is done then it should be documented. @@ -990,137 +1079,17 @@ overflow (including division by zero) which will result from executing these functions with an ill-conditioned matrix and thus raise Constraint_Error is sufficient. - -Implementations should not perform operations on mixed complex and real operands -by first converting the real operand to complex. See G.1.1(56,57). - - -G.3.3 Eigensystems of Matrices - -Eigenvalues and eigenvectors of a real symmetric matrix are computed using -subprograms in the child package Ada.Numerics.Generic_Real_Arrays.Eigen -whose declaration is: - -generic -package Ada.Numerics.Generic_Real_Arrays.Eigen is - pragma Pure(Eigen); - - -- Eigensystem of a real symmetric matrix - - 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 Argument_Error is -raised. - -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 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); - -end Ada.Numerics.Generic_Complex_Arrays.Eigen; - -The library package Numerics.Complex_Arrays.Eigen is declared pure and defines -the same subprograms as Numerics.Generic_Complex_Arrays.Eigen, except that -the Real_Vector and Real_Matrix types exported by Numerics.Real_Arrays are -systematically substituted for Real_Vector and Real_Matrix throughout. -Nongeneric equivalents for each of the other predefined floating point types -are defined similarly, with the names -Numerics.Short_Complex_Arrays.Eigen, Numerics.Long_Complex_Arrays.Eigen, etc. - -The effect of the subprograms is as described below. - -function Values(A : 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). - -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. - -For both subprograms, if the matrix A is not Hermitian then Argument_Error is -raised. - -Implementation Requirements - -The accuracy of these subprograms is implementation defined. - -Implementation Permissions - -The nongeneric equivalent packages may, but need not, be actual -instantiations of the generic package for the appropriate predefined type. - -Implementation Advice - -The test that a matrix is symmetric may be performed by using the equality -operator to compare the relevant components. The test that a matrix is -Hermitian may similarly use the equality operator to compare the real -components and negation followed by equality to compare the imaginary +The test that a matrix is Hermitian may use the equality operator to compare +the real components and negation followed by equality to compare the imaginary components (see G.2.1). -?? I hope that is good enough. I am not sure whether the RM implies that -negation is exact. If it doesn't then the paragraph above needs changing or -G.2.1 needs changing. +Implementations should not perform operations on mixed complex and real operands +by first converting the real operand to complex. See G.1.1(56, 57). !example +The packages are self-explanatory and so no example is provided. !discussion @@ -1139,12 +1108,16 @@ approach that we should add some basic facilities that are commonly required, not completely trivial to implement but on the other hand are mathematically well understood. + +The overall goal is thus twofold -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). +* to provide commonly required facilities for the user who is not a numerical + professional, +* 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 @@ -1174,13 +1147,20 @@ 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 +special operations to improve the speed and/or accuracy. We have specified +reasonable requirements in the exact mode, which are met by the canonical +implementation using a loop statement. This is because on some +hardware, built-in instructions which are fast actually lose accuracy. Note that the Fortran language standard recognizes the existence of inner product as a special case but imposes no accuracy requirements at all. +These requirements for inner product are based on the analysis in Chapter 3 +of Accuracy and Stability of Numerical Algorithms by N J Higham. A factor of +nearly 2 has been added as a precaution. In the complex case the factor is +increased by a factor of sqrt(2) in the mixed complex and real case to take +account of the hypotenuse effect and by a further factor of 2 in the pure +complex case since the number of terms is effectively doubled. + 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 @@ -1208,36 +1188,34 @@ DX := Solve(A, D); -and then correct X +and then correct X thus 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. +Implementations are recommended to do this automatically; it requires little +extra computation if LU decomposition is used internally. If they do such +refinement then it should be documented. + +A function for computing the determinant has been added since it is of some +help in deciding whether an array is ill-conditioned and therefore the results +of inversion might be suspect. Determinants are also useful in some statistical +calculations. The evaluation of determinants is very liable to overflow and many +such routines return a scaling power of 10 in order to keep the basic result +within range. For simplicity, it was decided not to do this since it is less +likely with matrices of low order; the user can of course scale the components +of the matrix if necessary. Similar functions have also been added for complex arrays. However, it was not deemed necessary to provide for mixed real and complex operands for Solve. In addition, subprograms have been added for the computation of eigenvalues -and vectors of real symmetric matrices and Hermitian matrices. Because these -are somewhat more specialized than the elementary operations they have been -placed in generic child packages such as Numerics.Generic_Real_Arrays.Eigen. -The individual subprograms are Values and Values_And_Vectors which means that -they can be invoked by statements such as - -V := Eigen.Values(A); - -There is no separate function for eigenvectors since it is unlikely that these -would be required without the eigenvalues. +and vectors of real symmetric matrices and Hermitian matrices. The subprograms +are Eigenvalues and Eigensystem. There is no separate subprogram for +eigenvectors since it is unlikely that these would be required without the +eigenvalues. The most common application is with real symmetric matrices and these have real -eigenvalues and vectors. The subprograms for these are in the real child -package. Applications include: moments of inertia, elasticity, +eigenvalues and vectors. Applications include: moments of inertia, elasticity, confidence regions and so on. A slight problem arises in deciding who should check that a matrix is symmetric, @@ -1261,12 +1239,16 @@ Note that errors such as when bounds of arrays do not match raise Constraint_Error by analogy with built-in array operations. +A third approach considered was for the user to supply a parameter indicating +which half of the matrix should be used to define it (the upper or lower +triangle). This avoids the need for any testing but it was considered bad +practice to be able to pass junk in the other half of the matrix. + The complex equivalent of real symmetric matrices are Hermitian matrices. Hermitian matrices are such that their transpose (that is with rows and columns interchanged) equals their complex conjugate (that is with the sign of imaginary -parts reversed). Hermitian matrices also have real eigenvalues and vectors. The -subprograms for these are in the complex child package. Applications include -Quantum Mechanics. +parts reversed). Hermitian matrices also have real eigenvalues. +Applications include Quantum Mechanics. Again we have placed the onus on the user to ensure that the matrix is Hermitian. The check in the package can then use strict equality for the real @@ -1300,7 +1282,7 @@ 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 +equations and determinant evaluation are to be performed using the same matrix. However, modern hardware is so fast that this would only seem to be necessary for very large sets of equations and these are not the target of these simple facilities. Moreover, adding explicit LU decomposition @@ -1314,6 +1296,1034 @@ not trivial for the user to program. These operations can of course be implemented as a binding to an implementation of part of the BLAS. +!corrigendum G.3(01) + +@dinsc + +Types and operations for the manipulation of real vectors and matrices are +provided in Generic_Real_Arrays, which is defined in G.3.1. Types and +operations for the manipulation of complex vectors and matrices are provided in +Generic_Complex_Arrays, which is defined in G.3.2. Both of these library units +are generic children of the predefined package Numerics (see A.5). Nongeneric +equivalents of these packages for each of the predefined floating point types +are also provided as children of Numerics. + +!corrigendum G.3.1(01) + +@dinsc + +@i<@s8<Static Semantics>> + +The generic library package Numerics.Generic_Real_Arrays has the following +declaration: + +@xcode<@b<generic> + @b<type> Real @b<is digits> <@>; +package Ada.Numerics.Generic_Real_Arrays is + @b<pragma> Pure(Generic_Real_Arrays); + + -- @ft<@i<Types>> + + @b<type> Real_Vector @b<is array> (Integer @b<range> <@>) @b<of> Real'Base; + @b<type> Real_Matrix @b<is array> (Integer @b<range> <@>, Integer @b<range> <@>) @b<of> Real'Base; + + -- @ft<@i<Subprograms for Real_Vector types>> + + -- @ft<@i<Real_Vector arithmetic operations>> + + @b<function> "+" (Right : Real_Vector) @b<return> Real_Vector; + @b<function> "-" (Right : Real_Vector) @b<return> Real_Vector; + @b<function> "@b<abs>" (Right : Real_Vector) @b<return> Real_Vector; + + @b<function> "+" (Left, Right : Real_Vector) @b<return> Real_Vector; + @b<function> "-" (Left, Right : Real_Vector) @b<return> Real_Vector; + + @b<function> "*" (Left, Right : Real_Vector) @b<return> Real'Base; + + -- @ft<@i<Real_Vector scaling operations>> + + @b<function> "*" (Left : Real'Base; Right : Real_Vector) @b<return> Real_Vector; + @b<function> "*" (Left : Real_Vector; Right : Real'Base) @b<return> Real_Vector; + @b<function> "/" (Left : Real_Vector; Right : Real'Base) @b<return> Real_Vector; + + -- @ft<@i<Other Real_Vector operations>> + + @b<function> Unit_Vector (Index : Integer; + Order : Positive; + First : Integer := 1) @b<return> Real_Vector; + + -- @ft<@i<Subprograms for Real_Matrix types>> + + -- @ft<@i<Real_Matrix arithmetic operations>> + + @b<function> "+" (Right : Real_Matrix) @b<return> Real_Matrix; + @b<function> "-" (Right : Real_Matrix) @b<return> Real_Matrix; + @b<function> "abs" (Right : Real_Matrix) @b<return> Real_Matrix; + @b<function> Transpose (X : Real_Matrix) @b<return> Real_Matrix; + + @b<function> "+" (Left, Right : Real_Matrix) @b<return> Real_Matrix; + @b<function> "-" (Left, Right : Real_Matrix) @b<return> Real_Matrix; + @b<function> "*" (Left, Right : Real_Matrix) @b<return> Real_Matrix; + + @b<function> "*" (Left, Right : Real_Vector) @b<return> Real_Matrix; + + @b<function> "*" (Left : Real_Vector; Right : Real_Matrix) @b<return> Real_Vector; + @b<function> "*" (Left : Real_Matrix; Right : Real_Vector) @b<return> Real_Vector; + + -- @ft<@i<Real_Matrix scaling operations>> + + @b<function> "*" (Left : Real'Base; Right : Real_Matrix) @b<return> Real_Matrix; + @b<function> "*" (Left : Real_Matrix; Right : Real'Base) @b<return> Real_Matrix; + @b<function> "/" (Left : Real_Matrix; Right : Real'Base) @b<return> Real_Matrix; + + -- @ft<@i<Real_Matrix inversion and related operations>> + + @b<function> Solve (A : Real_Matrix; X: Real_Vector) @b<return> Real_Vector; + @b<function> Solve (A, X : Real_Matrix) @b<return> Real_Matrix; + @b<function> Inverse (A : Real_Matrix) @b<return> Real_Matrix; + @b<function> Determinant (A : Real_Matrix) @b<return> Real'Base; + + -- @ft<@i<Eigenvalues and vectors of a real symmetric matrix>> + + @b<function> Eigenvalues(A : Real_Matrix) @b<return> Real_Vector; + + @b<procedure> Eigensystem(A : @b<in> Real_Matrix; + Values : @b<out> Real_Vector; + Vectors : @b<out> Real_Matrix); + + -- @ft<@i<Other Real_Matrix operations>> + + @b<function> Unit_Matrix (Order : Positive; + First_1, First_2 : Integer := 1) + @b<return> Real_Matrix; + +@b<private> + + -- @ft<@i<implementation-defined>> + +@b<end> Ada.Numerics.Generic_Real_Arrays;> + +The library package Numerics.Real_Arrays is declared pure and defines the +same types and subprograms as Numerics.Generic_Real_Arrays, except that +the predefined type Float is systematically substituted for Real'Base +throughout. Nongeneric equivalents for each of the other predefined floating +point types are defined similarly, with the names Numerics.Short_Real_Arrays, +Numerics.Long_Real_Arrays, etc. + +Two types are defined and exported by Ada.Numerics.Generic_Real_Arrays. The +composite type Real_Vector is provided to represent a vector with components of +type Real; it is defined as an unconstrained, one-dimensional array with an +index of type Integer. The composite type Real_Matrix is provided to represent +a matrix with components of type Real; it is defined as an unconstrained, +two-dimensional array with indices of type Integer. + +The effect of the various functions is as described below. In most cases the +functions are described in terms of corresponding scalar operations of the type +Real; any exception raised by those operations is propagated by the array +operation. Moreover the accuracy of the result for each individual component is +as defined for the scalar operation unless stated otherwise. + +In the case of those operations which are defined to involve an inner product, +Constraint_Error may be raised if an intermediate result is outside the range +of Real'Base even though the mathematical final result would not be. + +@xcode<@b<function> "+" (Right : Real_Vector) @b<return> Real_Vector; +@b<function> "-" (Right : Real_Vector) @b<return> Real_Vector; +@b<function> "@b<abs>" (Right : Real_Vector) @b<return> Real_Vector;> + +@xindent<Each operation returns the result of applying the corresponding operation of +the type Real to each component of Right. The index range of the result is +Right'Range.> + +@xcode<@b<function> "+" (Left, Right : Real_Vector) @b<return> Real_Vector; +@b<function> "-" (Left, Right : Real_Vector) @b<return> Real_Vector;> + +@xindent<Each operation returns the result of applying the corresponding operation of +the type Real to each component of Left and the matching component of Right. +The index range of the result is Left'Range. The exception Constraint_Error +is raised if Left'Length is not equal to Right'Length.> + +@xcode<@b<function> "*" (Left, Right : Real_Vector) @b<return> Real'Base;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real'Base; Right : Real_Vector) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Vector; Right : Real'Base) @b<return> Real_Vector; +@b<function> "/" (Left : Real_Vector; Right : Real'Base) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> Unit_Vector (Index : Integer; + Order : Positive; + First : Integer := 1) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> "+" (Right : Real_Matrix) @b<return> Real_Matrix; +@b<function> "-" (Right : Real_Matrix) @b<return> Real_Matrix; +@b<function> "abs" (Right : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<Each operation returns the result of applying the corresponding operation of +the type Real to each component of Right. The index ranges of the result are +those of Right.> + +@xcode<@b<function> Transpose (X : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "+" (Left, Right : Real_Matrix) @b<return> Real_Matrix; +@b<function> "-" (Left, Right : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<Each operation returns the result of applying the corresponding operation of +type Real to each component of Left and the matching component of Right. +The index ranges of the result are those of Left. The exception +Constraint_Error is raised if Left'Length(1) is not equal to Right'Length(1) or +Left'Length(2) is not equal to Right'Length(2).> + +@xcode<@b<function> "*" (Left, Right : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left, Right : Real_Vector) @b<return> Real_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Vector; Right : Real_Matrix) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Matrix; Right : Real_Vector) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real'Base; Right : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<This operation returns the result of multiplying each component of Right by the +scalar Left using the "*" operation of the type Real. The index ranges of the +matrix result are those of Right.> + +@xcode<@b<function> "*" (Left : Real_Matrix; Right : Real'Base) @b<return> Real_Matrix; +@b<function> "/" (Left : Real_Matrix; Right : Real'Base) @b<return> Real_Matrix;> + +@xindent<Each operation returns the result of applying the corresponding operation of +the type Real to each component of Left and to the scalar Right. The index +ranges of the matrix result are those of Left.> + +@xcode<@b<function> Solve (A : Real_Matrix; X: Real_Vector) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> Solve (A, X : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Inverse (A : Real_Matrix) @b<return> Real_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Determinant (A : Real_Matrix) @b<return> Real'Base;> + +@xindent<This function returns the determinant of the matrix A. Constraint_Error is +raised if A'Length(1) is not equal to A'Length(2).> + +@xcode<@b<function> Eigenvalues(A : Real_Matrix) @b<return> Real_Vector;> + +@xindent<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). The exception Argument_Error is raised if the matrix A +is not symmetric.> + +@xcode<@b<procedure> Eigensystem(A : @b<in> Real_Matrix; + Values : @b<out> Real_Vector; + Vectors : @b<out> Real_Matrix);> + +@xindent<This procedure computes both the eigenvalues and eigenvectors of the symmetric +matrix A. The out parameter Values is the same as that obtained by calling the +function Eigenvalues. The out parameter Vectors is a matrix whose columns are +the eigenvectors of the matrix A. The order of the columns corresponds to the +order of the eigenvalues. The eigenvectors are normalized and mutually +orthogonal (they are orthonormal), including when there are repeated +eigenvalues. 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 A is not symmetric.> + +@xcode<@b<function> Unit_Matrix (Order : Positive; + First_1, First_2 : Integer := 1) @b<return> Real_Matrix;> + +@xindent<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.> + +@i<@s8<Implementation Requirements>> + +Accuracy requirements for the subprograms Solve, Inverse, Determinant, +Eigenvalues and Eigensystem are implementation defined. + +For operations not involving an inner product, the accuracy +requirements are those of the corresponding operations of the type Real in +both the strict mode and the relaxed mode (see G.2). + +For operations involving an inner product, no requirements are specified in +the relaxed mode. In the strict mode the modulus of the absolute error of +the inner product X*Y shall not exceed g*abs(X)*abs(Y) where g is defined as +@xindent<g = X'Length * Machine_Radix**(1-Machine_Mantissa)> + +Implementations shall document any techniques used to reduce cancellation +errors such as extended precision arithmetic. + +@i<@s8<Implementation Permissions>> + +The nongeneric equivalent packages may, but need not, be actual +instantiations of the generic package for the appropriate predefined type. + +@i<@s8<Implementation Advice>> + +Implementations should implement the Solve and Inverse functions using +established techniques such as LU decomposition with row interchanges followed +by back and forward substitution. Implementations are recommended to refine the +result by performing an iteration on the residuals; if this is done then it +should be documented. + +It is not the intention that any special provision should be made to +determine whether a matrix is ill-conditioned or not. The naturally occurring +overflow (including division by zero) which will result from executing these +functions with an ill-conditioned matrix and thus raise Constraint_Error is +sufficient. + +The test that a matrix is symmetric may be performed by using the equality +operator to compare the relevant components. + +!corrigendum G.3.2(01) + +@dinsc + +@i<@s8<Static Semantics>> + +The generic library package Numerics.Generic_Complex_Arrays has the following +declaration: + +@xcode<@b<with> Ada.Numerics.Generic_Real_Arrays, Ada.Numerics.Generic_Complex_Types; +@b<generic> + @b<with package> Real_Arrays @b<is new> Ada.Numerics.Generic_Real_Arrays (<@>); + @b<use> Real_Arrays; + @b<with package> Complex_Types @b<is new> Ada.Numerics.Generic_Complex_Types (Real); + @b<use> Complex_Types; +@b<package> Ada.Numerics.Generic_Complex_Arrays @b<is> + @b<pragma> Pure(Generic_Complex_Arrays); + + -- @ft<@i<Types>> + + @b<type> Complex_Vector @b<is array> (Integer @b<range> <@>) @b<of> Complex; + @b<type> Complex_Matrix @b<is array> (Integer @b<range> <@>, + Integer @b<range> <@>) @b<of> Complex; + + -- @ft<@i<Subprograms for Complex_Vector types>> + + -- @ft<@i<Complex_Vector selection, conversion and composition operations>> + + @b<function> Re (X : Complex_Vector) @b<return> Real_Vector; + @b<function> Im (X : Complex_Vector) @b<return> Real_Vector; + + @b<procedure> Set_Re (X : @b<in out> Complex_Vector; + Re : @b<in> Real_Vector); + @b<procedure> Set_Im (X : @b<in out> Complex_Vector; + Im : @b<in> Real_Vector); + + @b<function> Compose_From_Cartesian (Re : Real_Vector) @b<return> Complex_Vector; + @b<function> Compose_From_Cartesian (Re, Im : Real_Vector) @b<return> Complex_Vector; + + @b<function> Modulus (X : Complex_Vector) @b<return> Real_Vector; + @b<function> "@b<abs>" (Right : Complex_Vector) @b<return> Real_Vector + @b<renames> Modulus; + @b<function> Argument (X : Complex_Vector) @b<return> Real_Vector; + @b<function> Argument (X : Complex_Vector; + Cycle : Real'Base) @b<return> Real_Vector; + + @b<function> Compose_From_Polar (Modulus, Argument : Real_Vector) + @b<return> Complex_Vector; + @b<function> Compose_From_Polar (Modulus, Argument : Real_Vector; + Cycle : Real'Base) + @b<return> Complex_Vector; + + -- @ft<@i<Complex_Vector arithmetic operations>> + + @b<function> "+" (Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "-" (Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> Conjugate (X : Complex_Vector) @b<return> Complex_Vector; + + @b<function> "+" (Left, Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "-" (Left, Right : Complex_Vector) @b<return> Complex_Vector; + + @b<function> "*" (Left, Right : Complex_Vector) @b<return> Complex; + + -- @ft<@i<Mixed Real_Vector and Complex_Vector arithmetic operations>> + + @b<function> "+" (Left : Real_Vector; + Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "+" (Left : Complex_Vector; + Right : Real_Vector) @b<return> Complex_Vector; + @b<function> "-" (Left : Real_Vector; + Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "-" (Left : Complex_Vector; + Right : Real_Vector) @b<return> Complex_Vector; + + @b<function> "*" (Left : Real_Vector; Right : Complex_Vector) @b<return> Complex; + @b<function> "*" (Left : Complex_Vector; Right : Real_Vector) @b<return> Complex; + + -- @ft<@i<Complex_Vector scaling operations>> + + @b<function> "*" (Left : Complex; + Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "*" (Left : Complex_Vector; + Right : Complex) @b<return> Complex_Vector; + @b<function> "/" (Left : Complex_Vector; + Right : Complex) @b<return> Complex_Vector; + + @b<function> "*" (Left : Real'Base; + Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "*" (Left : Complex_Vector; + Right : Real'Base) @b<return> Complex_Vector; + @b<function> "/" (Left : Complex_Vector; + Right : Real'Base) @b<return> Complex_Vector; + + -- @ft<@i<Other Complex_Vector operations>> + + @b<function> Unit_Vector (Index : Integer; + Order : Positive; + First : Integer := 1) @b<return> Complex_Vector; + + -- @ft<@i<Subprograms for Complex_Matrix types>> + + -- @ft<@i<Complex_Matrix selection, conversion and composition operations>> + + @b<function> Re (X : Complex_Matrix) @b<return> Real_Matrix; + @b<function> Im (X : Complex_Matrix) @b<return> Real_Matrix; + + @b<procedure> Set_Re (X : @b<in out> Complex_Matrix; + Re : @b<in> Real_Matrix); + @b<procedure> Set_Im (X : @b<in out> Complex_Matrix; + Im : @b<in> Real_Matrix); + + @b<function> Compose_From_Cartesian (Re : Real_Matrix) @b<return> Complex_Matrix; + @b<function> Compose_From_Cartesian (Re, Im : Real_Matrix) @b<return> Complex_Matrix; + + @b<function> Modulus (X : Complex_Matrix) @b<return> Real_Matrix; + @b<function> "@b<abs>" (Right : Complex_Matrix) @b<return> Real_Matrix + @b<renames> Modulus; + + @b<function> Argument (X : Complex_Matrix) @b<return> Real_Matrix; + @b<function> Argument (X : Complex_Matrix; + Cycle : Real'Base) @b<return> Real_Matrix; + + @b<function> Compose_From_Polar (Modulus, Argument : Real_Matrix) + @b<return> Complex_Matrix; + @b<function> Compose_From_Polar (Modulus, Argument : Real_Matrix; + Cycle : Real'Base) + @b<return> Complex_Matrix; + + -- @ft<@i<Complex_Matrix arithmetic operations>> + + @b<function> "+" (Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "-" (Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> Conjugate (X : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> Transpose (X : Complex_Matrix) @b<return> Complex_Matrix; + + @b<function> "+" (Left, Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "-" (Left, Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "*" (Left, Right : Complex_Matrix) @b<return> Complex_Matrix; + + @b<function> "*" (Left, Right : Complex_Vector) @b<return> Complex_Matrix; + + @b<function> "*" (Left : Complex_Vector; + Right : Complex_Matrix) @b<return> Complex_Vector; + @b<function> "*" (Left : Complex_Matrix; + Right : Complex_Vector) @b<return> Complex_Vector; + + -- @ft<@i<Mixed Real_Matrix and Complex_Matrix arithmetic operations>> + + @b<function> "+" (Left : Real_Matrix; + Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "+" (Left : Complex_Matrix; + Right : Real_Matrix) @b<return> Complex_Matrix; + @b<function> "-" (Left : Real_Matrix; + Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "-" (Left : Complex_Matrix; + Right : Real_Matrix) @b<return> Complex_Matrix; + @b<function> "*" (Left : Real_Matrix; + Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "*" (Left : Complex_Matrix; + Right : Real_Matrix) @b<return> Complex_Matrix; + + @b<function> "*" (Left : Real_Vector; + Right : Complex_Vector) @b<return> Complex_Matrix; + @b<function> "*" (Left : Complex_Vector; + Right : Real_Vector) @b<return> Complex_Matrix; + + @b<function> "*" (Left : Real_Vector; + Right : Complex_Matrix) @b<return> Complex_Vector; + @b<function> "*" (Left : Complex_Vector; + Right : Real_Matrix) @b<return> Complex_Vector; + @b<function> "*" (Left : Real_Matrix; + Right : Complex_Vector) @b<return> Complex_Vector; + @b<function> "*" (Left : Complex_Matrix; + Right : Real_Vector) @b<return> Complex_Vector; + + -- @ft<@i<Complex_Matrix scaling operations>> + + @b<function> "*" (Left : Complex; + Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "*" (Left : Complex_Matrix; + Right : Complex) @b<return> Complex_Matrix; + @b<function> "/" (Left : Complex_Matrix; + Right : Complex) @b<return> Complex_Matrix; + + @b<function> "*" (Left : Real'Base; + Right : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> "*" (Left : Complex_Matrix; + Right : Real'Base) @b<return> Complex_Matrix; + @b<function> "/" (Left : Complex_Matrix; + Right : Real'Base) @b<return> Complex_Matrix; + + -- @ft<@i<Complex_Matrix inversion and related operations>> + + @b<function> Solve (A : Complex_Matrix; X: Complex_Vector) @b<return> Complex_Vector; + @b<function> Solve (A, X : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> Inverse (A : Complex_Matrix) @b<return> Complex_Matrix; + @b<function> Determinant (A : Complex_Matrix) @b<return> Complex; + + -- @ft<@i<Eigenvalues and vectors of a Hermitian matrix>> + + @b<function> Eigenvalues(A : Complex_Matrix) @b<return> Real_Vector; + + @b<procedure> Eigensystem(A : @b<in> Complex_Matrix; + Values : @b<out> Real_Vector; + Vectors : @b<out> Complex_Matrix); + + -- @ft<@i<Other Complex_Matrix operations>> + + @b<function> Unit_Matrix (Order : Positive; + First_1, First_2 : Integer := 1) + @b<return> Complex_Matrix; + +@b<end> Ada.Numerics.Generic_Complex_Arrays;> + +The library package Numerics.Complex_Arrays is declared pure and defines +the same types and subprograms as Numerics.Generic_Complex_Arrays, except +that the predefined type Float is systematically substituted for Real'Base, +and the Real_Vector and Real_Matrix types exported by Numerics.Real_Arrays +are systematically substituted for Real_Vector and Real_Matrix, and the +Complex type exported by Numerics.Complex_Types is systematically +substituted for Complex, throughout. Nongeneric equivalents for each of +the other predefined floating point types are defined similarly, with the +names Numerics.Short_Complex_Arrays, Numerics.Long_Complex_Arrays, etc. + +Two types are defined and exported by Ada.Numerics.Generic_Complex_Arrays. +The composite type Complex_Vector is provided to represent a vector with +components of type Complex; it is defined as an unconstrained +one-dimensional array with an index of type Integer. The composite type +Complex_Matrix is provided to represent a matrix with components of type +Complex; it is defined as an unconstrained, two-dimensional array with +indices of type Integer. + +The effect of the various subprograms is as described below. In many cases +they are described in terms of corresponding scalar operations in +Numerics.Generic_Complex_Types. Any exception raised by those operations is +propagated by the array subprogram. Moreover any constraints on the parameters +and the accuracy of the result for each individual component 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. + +@xcode<@b<function> Re (X : Complex_Vector) @b<return> Real_Vector; +@b<function> Im (X : Complex_Vector) @b<return> Real_Vector;> + +@xindent<Each function returns a vector of the specified cartesian components of X. +The index range of the result is X'Range.> + +@xcode<@b<procedure> Set_Re (X : @b<in out> Complex_Vector; Re : @b<in> Real_Vector); +@b<procedure> Set_Im (X : @b<in out> Complex_Vector; Im : @b<in> Real_Vector);> + +@xindent<Each procedure replaces the specified (cartesian) component of each of +the components of X by the value of the matching component of Re or Im; +the other (cartesian) component of each of the components is unchanged. +The exception Constraint_Error is raised if X'Length is not equal to +Re'Length or Im'Length.> + +@xcode<@b<function> Compose_From_Cartesian (Re : Real_Vector) @b<return> Complex_Vector; +@b<function> Compose_From_Cartesian (Re, Im : Real_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> Modulus (X : Complex_Vector) @b<return> Real_Vector; +@b<function> "abs" (Right : Complex_Vector) @b<return> Real_Vector @b<renames> Modulus; +@b<function> Argument (X : Complex_Vector) @b<return> Real_Vector; +@b<function> Argument (X : Complex_Vector; + Cycle : Real'Base) @b<return> Real_Vector;> + +@xindent<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.> + +@xcode<@b<function> Compose_From_Polar (Modulus, Argument : Real_Vector) + @b<return> Complex_Vector; +@b<function> Compose_From_Polar (Modulus, Argument : Real_Vector; Cycle : Real'Base) + @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "+" (Right : Complex_Vector) @b<return> Complex_Vector; +@b<function> "-" (Right : Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> Conjugate (X : Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "+" (Left, Right : Complex_Vector) @b<return> Complex_Vector; +@b<function> "-" (Left, Right : Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left, Right : Complex_Vector) @b<return> Complex;> + +@xindent<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.> + +@xcode<@b<function> "+" (Left : Real_Vector; + Right : Complex_Vector) @b<return> Complex_Vector; +@b<function> "+" (Left : Complex_Vector; + Right : Real_Vector) @b<return> Complex_Vector; +@b<function> "-" (Left : Real_Vector; + Right : Complex_Vector) @b<return> Complex_Vector; +@b<function> "-" (Left : Complex_Vector; + Right : Real_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Vector; Right : Complex_Vector) @b<return> Complex; +@b<function> "*" (Left : Complex_Vector; Right : Real_Vector) @b<return> Complex;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex; Right : Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex_Vector; Right : Complex) @b<return> Complex_Vector; +@b<function> "/" (Left : Complex_Vector; Right : Complex) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real'Base; Right : Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex_Vector; Right : Real'Base) @b<return> Complex_Vector; +@b<function> "/" (Left : Complex_Vector; Right : Real'Base) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> Unit_Vector (Index : Integer; + Order : Positive; + First : Integer := 1) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> Re (X : Complex_Matrix) @b<return> Real_Matrix; +@b<function> Im (X : Complex_Matrix) @b<return> Real_Matrix;> + +@xindent<Each function returns a matrix of the specified cartesian components of X. +The index ranges of the result are those of X.> + +@xcode<@b<procedure> Set_Re (X : @b<in out> Complex_Matrix; Re : @b<in> Real_Matrix); +@b<procedure> Set_Im (X : @b<in out> Complex_Matrix; Im : @b<in> Real_Matrix);> + +@xindent<Each procedure replaces the specified (cartesian) component of each of the +components of X by the value of the matching component of Re or Im; the other +(cartesian) component of each of the components is unchanged. The exception +Constraint_Error is raised if X'Length(1) is not equal to Re'Length(1) or +Im'Length(1) or if X'Length(2) is not equal to Re'Length(2) or Im'Length(2).> + +@xcode<@b<function> Compose_From_Cartesian (Re : Real_Matrix) @b<return> Complex_Matrix; +@b<function> Compose_From_Cartesian (Re, Im : Real_Matrix) @b<return> Complex_Matrix;> + +@xindent<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).> + +@xcode<@b<function> Modulus (X : Complex_Matrix) @b<return> Real_Matrix; +@b<function> "@b<abs>" (Right : Complex_Matrix) @b<return> Real_Matrix @b<renames> Modulus; +@b<function> Argument (X : Complex_Matrix) @b<return> Real_Matrix; +@b<function> Argument (X : Complex_Matrix; + Cycle : Real'Base) @b<return> Real_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Compose_From_Polar (Modulus, Argument : Real_Matrix) + @b<return> Complex_Matrix; +@b<function> Compose_From_Polar (Modulus, Argument : Real_Matrix; + Cycle : Real'Base) + @b<return> Complex_Matrix;> + +@xindent<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).> + +@xcode<@b<function> "+" (Right : Complex_Matrix) @b<return> Complex_Matrix; +@b<function> "-" (Right : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Conjugate (X : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Transpose (X : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "+" (Left, Right : Complex_Matrix) @b<return> Complex_Matrix; +@b<function> "-" (Left, Right : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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).> + +@xcode<@b<function> "*" (Left, Right : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left, Right : Complex_Vector) @b<return> Complex_Matrix;> + +@xindent<This operation returns the outer product of a (column) vector Left by a (row) +vector Right using the appropriate operation "*" in +Numerics.Generic_Complex_Types for computing the individual components. +The first and second index ranges of the matrix result are Left'Range and +Right'Range respectively.> + +@xcode<@b<function> "*" (Left : Complex_Vector; + Right : Complex_Matrix) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex_Matrix; + Right : Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "+" (Left : Real_Matrix; + Right : Complex_Matrix) @b<return> Complex_Matrix; +@b<function> "+" (Left : Complex_Matrix; + Right : Real_Matrix) @b<return> Complex_Matrix; +@b<function> "-" (Left : Real_Matrix; + Right : Complex_Matrix) @b<return> Complex_Matrix; +@b<function> "-" (Left : Complex_Matrix; + Right : Real_Matrix) @b<return> Complex_Matrix;> + +@xindent<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).> + +@xcode<@b<function> "*" (Left : Real_Matrix; + Right : Complex_Matrix) @b<return> Complex_Matrix; +@b<function> "*" (Left : Complex_Matrix; + Right : Real_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Vector; + Right : Complex_Vector) @b<return> Complex_Matrix; +@b<function> "*" (Left : Complex_Vector; + Right : Real_Vector) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Vector; + Right : Complex_Matrix) @b<return> Complex_Vector; +@b<function> "*" (Left : Complex_Vector; + Right : Real_Matrix) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real_Matrix; + Right : Complex_Vector) @b<return> Complex_Vector; +@b<function> "*" (Left : Complex_Matrix; + Right : Real_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex; Right : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex_Matrix; Right : Complex) @b<return> Complex_Matrix; +@b<function> "/" (Left : Complex_Matrix; Right : Complex) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Real'Base; Right : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> "*" (Left : Complex_Matrix; Right : Real'Base) @b<return> Complex_Matrix; +@b<function> "/" (Left : Complex_Matrix; Right : Real'Base) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Solve (A : Complex_Matrix; X: Complex_Vector) @b<return> Complex_Vector;> + +@xindent<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.> + +@xcode<@b<function> Solve (A, X : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Inverse (A : Complex_Matrix) @b<return> Complex_Matrix;> + +@xindent<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.> + +@xcode<@b<function> Determinant (A : Complex_Matrix) @b<return> Complex;> + +@xindent<This function returns the determinant of the matrix A. Constraint_Error is +raised if A'Length(1) is not equal to A'Length(2).> + +@xcode<@b<function> Eigenvalues(A : Complex_Matrix) @b<return> Real_Vector;> + +@xindent<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 A +is not Hermitian.> + +@xcode<@b<procedure> Eigensystem(A : @b<in> Complex_Matrix; + Values : @b<out> Real_Vector; + Vectors : @b<out> Complex_Matrix);> + +@xindent<This procedure computes both the eigenvalues and eigenvectors of the Hermitian +matrix A. The out parameter Values is the same as that obtained by calling the +function Eigenvalues. The out parameter Vectors is a matrix whose columns are +the eigenvectors of the matrix A. The order of the columns corresponds to the +order of the eigenvalues. The eigenvectors are mutually orthonormal, +including when there are repeated eigenvalues. 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 A is not Hermitian.> + +@xcode<@b<function> Unit_Matrix (Order : Positive; + First_1, First_2 : Integer := 1) + @b<return> Complex_Matrix;> + +@xindent<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.> + +@i<@s8<Implementation Requirements>> + +Accuracy requirements for the subprograms Solve, Inverse, Determinant, +Eigenvalues and Eigensystem are implementation defined. + +For operations not involving an inner product, the accuracy requirements are +those of the corresponding operations of the type Real'Base and Complex in both +the strict mode and the relaxed mode (see G.2). + +For operations involving an inner product, no requirements are specified in +the relaxed mode. In the strict mode the modulus of the absolute error of the +inner product X*Y shall not exceed g*abs(X)*abs(Y) where g is defined as +@xindent<g = sqrt(2.0) * X'Length * Machine_Radix**(1-Machine_Mantissa) for mixed + complex and real operands> + +@xindent<g = sqrt(8.0) * X'Length * Machine_Radix**(1-Machine_Mantissa) for two + complex operands> + +Implementations shall document any techniques used to reduce cancellation +errors such as extended precision arithmetic. + +@i<@s8<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. + +@i<@s8<Implementation Advice>> + +Implementations should implement the Solve and Inverse functions using +established techniques. Implementations are recommended to refine the result by +performing an iteration on the residuals; if this is done then it should be +documented. + +It is not the intention that any special provision should be made to +determine whether a matrix is ill-conditioned or not. The naturally occurring +overflow (including division by zero) which will result from executing these +functions with an ill-conditioned matrix and thus raise Constraint_Error is +sufficient. + +The test that a matrix is Hermitian may use the equality operator to compare +the real components and negation followed by equality to compare the imaginary +components (see G.2.1). + +Implementations should not perform operations on mixed complex and real operands +by first converting the real operand to complex. See G.1.1. + + + + + !ACATS Test ACATS test(s) need to be created. @@ -1604,7 +2614,7 @@ would be easy. Curve and surface fitting. There is lots of choice here. -Cubic splines are fashionable. Chebychev polynolmials are a +Cubic splines are fashionable. Chebychev polynomials are a long established technique. But both require significant user knowledge and it is easy to misuse them.

Questions? Ask the ACAA Technical Agent