Next: , Previous: G.3, Up: G.3


G.3.1 Real Vectors and Matrices

Static Semantics

1/2
The generic library package Numerics.Generic_Real_Arrays has the following declaration:

2/2

     generic
        type Real is digits <>;
     package Ada.Numerics.Generic_Real_Arrays is   pragma Pure(Generic_Real_Arrays);

3/2

        −− Types

4/2

        type Real_Vector is array (Integer range <>) of Real'Base;
        type Real_Matrix is array (Integer range <>, Integer range <>)
                                                        of Real'Base;

5/2

        −− Subprograms for Real_Vector types

6/2

        −− Real_Vector arithmetic operations

7/2

        function "+"   (Right Real_Vector)       return Real_Vector;
        function "−"   (Right Real_Vector)       return Real_Vector;
        function "abs" (Right Real_Vector)       return Real_Vector;

8/2

        function "+"   (Left, Right Real_Vector) return Real_Vector;
        function "−"   (Left, Right Real_Vector) return Real_Vector;

9/2

        function "*"   (Left, Right Real_Vector) return Real'Base;

10/2

        function "abs" (Right Real_Vector)       return Real'Base;

11/2

        −− Real_Vector scaling operations

12/2

        function "*" (Left Real'Base;   Right Real_Vector)
           return Real_Vector;
        function "*" (Left Real_Vector; Right Real'Base)
           return Real_Vector;
        function "/" (Left Real_Vector; Right Real'Base)
           return Real_Vector;

13/2

        −− Other Real_Vector operations

14/2

        function Unit_Vector (Index Integer;
                              Order Positive;
                              First Integer := 1) return Real_Vector;

15/2

        −− Subprograms for Real_Matrix types

16/2

        −− Real_Matrix arithmetic operations

17/2

        function "+"       (Right Real_Matrix) return Real_Matrix;
        function "−"       (Right Real_Matrix) return Real_Matrix;
        function "abs"     (Right Real_Matrix) return Real_Matrix;
        function Transpose (X     Real_Matrix) return Real_Matrix;

18/2

        function "+" (Left, Right Real_Matrix) return Real_Matrix;
        function "−" (Left, Right Real_Matrix) return Real_Matrix;
        function "*" (Left, Right Real_Matrix) return Real_Matrix;

19/2

        function "*" (Left, Right Real_Vector) return Real_Matrix;

20/2

        function "*" (Left Real_Vector; Right Real_Matrix)
           return Real_Vector;
        function "*" (Left Real_Matrix; Right Real_Vector)
           return Real_Vector;

21/2

        −− Real_Matrix scaling operations

22/2

        function "*" (Left Real'Base;   Right Real_Matrix)
           return Real_Matrix;
        function "*" (Left Real_Matrix; Right Real'Base)
           return Real_Matrix;
        function "/" (Left Real_Matrix; Right Real'Base)
           return Real_Matrix;

23/2

        −− Real_Matrix inversion and related operations

24/2

        function Solve (A Real_Matrix; Real_Vector) return Real_Vector;
        function Solve (A, Real_Matrix) return Real_Matrix;
        function Inverse (A Real_Matrix) return Real_Matrix;
        function Determinant (A Real_Matrix) return Real'Base;

25/2

        −− Eigenvalues and vectors of real symmetric matrix

26/2

        function Eigenvalues (A Real_Matrix) return Real_Vector;

27/2

        procedure Eigensystem (A       in  Real_Matrix;
                               Values  out Real_Vector;
                               Vectors out Real_Matrix);

28/2

        −− Other Real_Matrix operations

29/2

        function Unit_Matrix (Order            Positive;
                              First_1, First_2 Integer := 1)
                                                 return Real_Matrix;

30/2

     end Ada.Numerics.Generic_Real_Arrays;

31/2
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.

32/2
Two types are defined and exported by 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.

33/2
The effect of the various subprograms is as described below. In most cases the subprograms 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.

34/2
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.35/2

     function "+"   (Right Real_Vector) return Real_Vector;
     function "−"   (Right Real_Vector) return Real_Vector;
     function "abs" (Right Real_Vector) return Real_Vector;

36/2

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.

37/2

     function "+" (Left, Right Real_Vector) return Real_Vector;
     function "−" (Left, Right Real_Vector) return Real_Vector;

38/2

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. Constraint_Error is raised if Left'Length is not equal to Right'Length.

39/2

     function "*" (Left, Right Real_Vector) return Real'Base;

40/2

This operation returns the inner product of Left and Right. Constraint_Error is raised if Left'Length is not equal to Right'Length. This operation involves an inner product.

41/2

     function "abs" (Right Real_Vector) return Real'Base;

42/2

This operation returns the L2−norm of Right (the square root of the inner product of the vector with itself).

43/2

     function "*" (Left Real'Base; Right Real_Vector) return Real_Vector;

44/2

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.

45/2

     function "*" (Left Real_Vector; Right Real'Base) return Real_Vector;
     function "/" (Left Real_Vector; Right Real'Base) return Real_Vector;

46/2

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.

47/2

     function Unit_Vector (Index Integer;
                           Order Positive;
                           First Integer := 1) return Real_Vector;

48/2

This function returns a unit vectorwith 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. Constraint_Error is raised if Index < First, Index > First + Order −− 1 or if First + Order −− 1 > Integer'Last.

49/2

     function "+"   (Right Real_Matrix) return Real_Matrix;
     function "−"   (Right Real_Matrix) return Real_Matrix;
     function "abs" (Right Real_Matrix) return Real_Matrix;

50/2

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.

51/2

     function Transpose (X Real_Matrix) return Real_Matrix;

52/2

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.

53/2

     function "+" (Left, Right Real_Matrix) return Real_Matrix;
     function "−" (Left, Right Real_Matrix) return Real_Matrix;

54/2

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 ranges of the result are those of Left. 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).

55/2

     function "*" (Left, Right Real_Matrix) return Real_Matrix;

56/2

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. Constraint_Error is raised if Left'Length(2) is not equal to Right'Length(1). This operation involves inner products.

57/2

     function "*" (Left, Right Real_Vector) return Real_Matrix;

58/2

This operation returns the outer product of a (column) vector Left by a (row) vector Right using the operation "*" of the type Real for computing the individual components. The first and second index ranges of the result are Left'Range and Right'Range respectively.

59/2

     function "*" (Left Real_Vector; Right Real_Matrix) return Real_Vector;

60/2

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). Constraint_Error is raised if Left'Length is not equal to Right'Length(1). This operation involves inner products.

61/2

     function "*" (Left Real_Matrix; Right Real_Vector) return Real_Vector;

62/2

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). Constraint_Error is raised if Left'Length(2) is not equal to Right'Length. This operation involves inner products.

63/2

     function "*" (Left Real'Base; Right Real_Matrix) return Real_Matrix;

64/2

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 result are those of Right.

65/2

     function "*" (Left Real_Matrix; Right Real'Base) return Real_Matrix;
     function "/" (Left Real_Matrix; Right Real'Base) return Real_Matrix;

66/2

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 result are those of Left.

67/2

     function Solve (A Real_Matrix; Real_Vector) return Real_Vector;

68/2

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 A'Range(2). 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.

69/2

     function Solve (A, Real_Matrix) return Real_Matrix;

70/2

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 A'Range(2) and X'Range(2). 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.

71/2

     function Inverse (A Real_Matrix) return Real_Matrix;

72/2

This function returns a matrix B such that A * B is (nearly) equal to the unit matrix. The index ranges of the result are A'Range(2) and A'Range(1). 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.

73/2

     function Determinant (A Real_Matrix) return Real'Base;

74/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).

75/2

     function Eigenvalues(A Real_Matrix) return Real_Vector;

76/2

This function returns the eigenvalues of the symmetric matrix A as a vector sorted into order with the largest first. 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). Argument_Error is raised if the matrix A is not symmetric.

77/2

     procedure Eigensystem(A       in  Real_Matrix;
                           Values  out Real_Vector;
                           Vectors out Real_Matrix);

78/2

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. 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. Argument_Error is raised if the matrix A is not symmetric.

79/2

     function Unit_Matrix (Order            Positive;
                           First_1, First_2 Integer := 1) return Real_Matrix;

80/2

This function returns a square unit matrixwith 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. Constraint_Error is raised if First_1 + Order −− 1 > Integer'Last or First_2 + Order −− 1 > Integer'Last.
Implementation Requirements

81/2
Accuracy requirements for the subprograms Solve, Inverse, Determinant, Eigenvalues and Eigensystem are implementation defined.

82/2
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).

83/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

84/2

X'Length Real'Machine_Radix**(1 −− Real'Model_Mantissa)

85/2
For the L2−norm, no accuracy requirements are specified in the relaxed mode. In the strict mode the relative error on the norm shall not exceed g / 2.0 + 3.0 * Real'Model_Epsilon where g is defined as above.

Documentation Requirements

86/2
Implementations shall document any techniques used to reduce cancellation errors such as extended precision arithmetic.

Implementation Permissions

87/2
The nongeneric equivalent packages may, but need not, be actual instantiations of the generic package for the appropriate predefined type.

Implementation Advice

88/2
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.

89/2
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.

90/2
The test that a matrix is symmetric should be performed by using the equality operator to compare the relevant components.