------------------------------------------------------------------------------ -- -- -- GNAT RUN-TIME COMPONENTS -- -- -- -- ADA.NUMERICS.GENERIC_REAL_ARRAYS -- -- -- -- B o d y -- -- -- -- Copyright (C) 2006-2019, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- -- ware Foundation; either version 3, or (at your option) any later ver- -- -- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- -- or FITNESS FOR A PARTICULAR PURPOSE. -- -- -- -- As a special exception under Section 7 of GPL version 3, you are granted -- -- additional permissions described in the GCC Runtime Library Exception, -- -- version 3.1, as published by the Free Software Foundation. -- -- -- -- You should have received a copy of the GNU General Public License and -- -- a copy of the GCC Runtime Library Exception along with this program; -- -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- -- . -- -- -- -- GNAT was originally developed by the GNAT team at New York University. -- -- Extensive contributions were provided by Ada Core Technologies Inc. -- -- -- ------------------------------------------------------------------------------ -- This version of Generic_Real_Arrays avoids the use of BLAS and LAPACK. One -- reason for this is new Ada 2012 requirements that prohibit algorithms such -- as Strassen's algorithm, which may be used by some BLAS implementations. In -- addition, some platforms lacked suitable compilers to compile the reference -- BLAS/LAPACK implementation. Finally, on some platforms there are more -- floating point types than supported by BLAS/LAPACK. with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers; with System; use System; with System.Generic_Array_Operations; use System.Generic_Array_Operations; package body Ada.Numerics.Generic_Real_Arrays is package Ops renames System.Generic_Array_Operations; function Is_Non_Zero (X : Real'Base) return Boolean is (X /= 0.0); procedure Back_Substitute is new Ops.Back_Substitute (Scalar => Real'Base, Matrix => Real_Matrix, Is_Non_Zero => Is_Non_Zero); function Diagonal is new Ops.Diagonal (Scalar => Real'Base, Vector => Real_Vector, Matrix => Real_Matrix); procedure Forward_Eliminate is new Ops.Forward_Eliminate (Scalar => Real'Base, Real => Real'Base, Matrix => Real_Matrix, Zero => 0.0, One => 1.0); procedure Swap_Column is new Ops.Swap_Column (Scalar => Real'Base, Matrix => Real_Matrix); procedure Transpose is new Ops.Transpose (Scalar => Real'Base, Matrix => Real_Matrix); function Is_Symmetric (A : Real_Matrix) return Boolean is (Transpose (A) = A); -- Return True iff A is symmetric, see RM G.3.1 (90). function Is_Tiny (Value, Compared_To : Real) return Boolean is (abs Compared_To + 100.0 * abs (Value) = abs Compared_To); -- Return True iff the Value is much smaller in magnitude than the least -- significant digit of Compared_To. procedure Jacobi (A : Real_Matrix; Values : out Real_Vector; Vectors : out Real_Matrix; Compute_Vectors : Boolean := True); -- Perform Jacobi's eigensystem algorithm on real symmetric matrix A function Length is new Square_Matrix_Length (Real'Base, Real_Matrix); -- Helper function that raises a Constraint_Error is the argument is -- not a square matrix, and otherwise returns its length. procedure Rotate (X, Y : in out Real; Sin, Tau : Real); -- Perform a Givens rotation procedure Sort_Eigensystem (Values : in out Real_Vector; Vectors : in out Real_Matrix); -- Sort Values and associated Vectors by decreasing absolute value procedure Swap (Left, Right : in out Real); -- Exchange Left and Right function Sqrt is new Ops.Sqrt (Real); -- Instant a generic square root implementation here, in order to avoid -- instantiating a complete copy of Generic_Elementary_Functions. -- Speed of the square root is not a big concern here. ------------ -- Rotate -- ------------ procedure Rotate (X, Y : in out Real; Sin, Tau : Real) is Old_X : constant Real := X; Old_Y : constant Real := Y; begin X := Old_X - Sin * (Old_Y + Old_X * Tau); Y := Old_Y + Sin * (Old_X - Old_Y * Tau); end Rotate; ---------- -- Swap -- ---------- procedure Swap (Left, Right : in out Real) is Temp : constant Real := Left; begin Left := Right; Right := Temp; end Swap; -- Instantiating the following subprograms directly would lead to -- name clashes, so use a local package. package Instantiations is function "+" is new Vector_Elementwise_Operation (X_Scalar => Real'Base, Result_Scalar => Real'Base, X_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "+"); function "+" is new Matrix_Elementwise_Operation (X_Scalar => Real'Base, Result_Scalar => Real'Base, X_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "+"); function "+" is new Vector_Vector_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Right_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "+"); function "+" is new Matrix_Matrix_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Matrix => Real_Matrix, Right_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "+"); function "-" is new Vector_Elementwise_Operation (X_Scalar => Real'Base, Result_Scalar => Real'Base, X_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "-"); function "-" is new Matrix_Elementwise_Operation (X_Scalar => Real'Base, Result_Scalar => Real'Base, X_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "-"); function "-" is new Vector_Vector_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Right_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "-"); function "-" is new Matrix_Matrix_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Matrix => Real_Matrix, Right_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "-"); function "*" is new Scalar_Vector_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Right_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "*"); function "*" is new Scalar_Matrix_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Right_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "*"); function "*" is new Vector_Scalar_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "*"); function "*" is new Matrix_Scalar_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "*"); function "*" is new Outer_Product (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Right_Vector => Real_Vector, Matrix => Real_Matrix); function "*" is new Inner_Product (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Right_Vector => Real_Vector, Zero => 0.0); function "*" is new Matrix_Vector_Product (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Matrix => Real_Matrix, Right_Vector => Real_Vector, Result_Vector => Real_Vector, Zero => 0.0); function "*" is new Vector_Matrix_Product (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Matrix => Real_Matrix, Result_Vector => Real_Vector, Zero => 0.0); function "*" is new Matrix_Matrix_Product (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Matrix => Real_Matrix, Right_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Zero => 0.0); function "/" is new Vector_Scalar_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "/"); function "/" is new Matrix_Scalar_Elementwise_Operation (Left_Scalar => Real'Base, Right_Scalar => Real'Base, Result_Scalar => Real'Base, Left_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "/"); function "abs" is new L2_Norm (X_Scalar => Real'Base, Result_Real => Real'Base, X_Vector => Real_Vector, "abs" => "+"); -- While the L2_Norm by definition uses the absolute values of the -- elements of X_Vector, for real values the subsequent squaring -- makes this unnecessary, so we substitute the "+" identity function -- instead. function "abs" is new Vector_Elementwise_Operation (X_Scalar => Real'Base, Result_Scalar => Real'Base, X_Vector => Real_Vector, Result_Vector => Real_Vector, Operation => "abs"); function "abs" is new Matrix_Elementwise_Operation (X_Scalar => Real'Base, Result_Scalar => Real'Base, X_Matrix => Real_Matrix, Result_Matrix => Real_Matrix, Operation => "abs"); function Solve is new Matrix_Vector_Solution (Real'Base, 0.0, Real_Vector, Real_Matrix); function Solve is new Matrix_Matrix_Solution (Real'Base, 0.0, Real_Matrix); function Unit_Matrix is new Generic_Array_Operations.Unit_Matrix (Scalar => Real'Base, Matrix => Real_Matrix, Zero => 0.0, One => 1.0); function Unit_Vector is new Generic_Array_Operations.Unit_Vector (Scalar => Real'Base, Vector => Real_Vector, Zero => 0.0, One => 1.0); end Instantiations; --------- -- "+" -- --------- function "+" (Right : Real_Vector) return Real_Vector renames Instantiations."+"; function "+" (Right : Real_Matrix) return Real_Matrix renames Instantiations."+"; function "+" (Left, Right : Real_Vector) return Real_Vector renames Instantiations."+"; function "+" (Left, Right : Real_Matrix) return Real_Matrix renames Instantiations."+"; --------- -- "-" -- --------- function "-" (Right : Real_Vector) return Real_Vector renames Instantiations."-"; function "-" (Right : Real_Matrix) return Real_Matrix renames Instantiations."-"; function "-" (Left, Right : Real_Vector) return Real_Vector renames Instantiations."-"; function "-" (Left, Right : Real_Matrix) return Real_Matrix renames Instantiations."-"; --------- -- "*" -- --------- -- Scalar multiplication function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector renames Instantiations."*"; function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector renames Instantiations."*"; function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix renames Instantiations."*"; function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix renames Instantiations."*"; -- Vector multiplication function "*" (Left, Right : Real_Vector) return Real'Base renames Instantiations."*"; function "*" (Left, Right : Real_Vector) return Real_Matrix renames Instantiations."*"; function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector renames Instantiations."*"; function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector renames Instantiations."*"; -- Matrix Multiplication function "*" (Left, Right : Real_Matrix) return Real_Matrix renames Instantiations."*"; --------- -- "/" -- --------- function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector renames Instantiations."/"; function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix renames Instantiations."/"; ----------- -- "abs" -- ----------- function "abs" (Right : Real_Vector) return Real'Base renames Instantiations."abs"; function "abs" (Right : Real_Vector) return Real_Vector renames Instantiations."abs"; function "abs" (Right : Real_Matrix) return Real_Matrix renames Instantiations."abs"; ----------------- -- Determinant -- ----------------- function Determinant (A : Real_Matrix) return Real'Base is M : Real_Matrix := A; B : Real_Matrix (A'Range (1), 1 .. 0); R : Real'Base; begin Forward_Eliminate (M, B, R); return R; end Determinant; ----------------- -- Eigensystem -- ----------------- procedure Eigensystem (A : Real_Matrix; Values : out Real_Vector; Vectors : out Real_Matrix) is begin Jacobi (A, Values, Vectors, Compute_Vectors => True); Sort_Eigensystem (Values, Vectors); end Eigensystem; ----------------- -- Eigenvalues -- ----------------- function Eigenvalues (A : Real_Matrix) return Real_Vector is begin return Values : Real_Vector (A'Range (1)) do declare Vectors : Real_Matrix (1 .. 0, 1 .. 0); begin Jacobi (A, Values, Vectors, Compute_Vectors => False); Sort_Eigensystem (Values, Vectors); end; end return; end Eigenvalues; ------------- -- Inverse -- ------------- function Inverse (A : Real_Matrix) return Real_Matrix is (Solve (A, Unit_Matrix (Length (A), First_1 => A'First (2), First_2 => A'First (1)))); ------------ -- Jacobi -- ------------ procedure Jacobi (A : Real_Matrix; Values : out Real_Vector; Vectors : out Real_Matrix; Compute_Vectors : Boolean := True) is -- This subprogram uses Carl Gustav Jacob Jacobi's iterative method -- for computing eigenvalues and eigenvectors and is based on -- Rutishauser's implementation. -- The given real symmetric matrix is transformed iteratively to -- diagonal form through a sequence of appropriately chosen elementary -- orthogonal transformations, called Jacobi rotations here. -- The Jacobi method produces a systematic decrease of the sum of the -- squares of off-diagonal elements. Convergence to zero is quadratic, -- both for this implementation, as for the classic method that doesn't -- use row-wise scanning for pivot selection. -- The numerical stability and accuracy of Jacobi's method make it the -- best choice here, even though for large matrices other methods will -- be significantly more efficient in both time and space. -- While the eigensystem computations are absolutely foolproof for all -- real symmetric matrices, in presence of invalid values, or similar -- exceptional situations it might not. In such cases the results cannot -- be trusted and Constraint_Error is raised. -- Note: this implementation needs temporary storage for 2 * N + N**2 -- values of type Real. Max_Iterations : constant := 50; N : constant Natural := Length (A); subtype Square_Matrix is Real_Matrix (1 .. N, 1 .. N); -- In order to annihilate the M (Row, Col) element, the -- rotation parameters Cos and Sin are computed as -- follows: -- Theta = Cot (2.0 * Phi) -- = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col)) -- Then Tan (Phi) as the smaller root (in modulus) of -- T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large) function Compute_Tan (Theta : Real) return Real is (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta)); function Compute_Tan (P, H : Real) return Real is (if Is_Tiny (P, Compared_To => H) then P / H else Compute_Tan (Theta => H / (2.0 * P))); function Sum_Strict_Upper (M : Square_Matrix) return Real; -- Return the sum of all elements in the strict upper triangle of M ---------------------- -- Sum_Strict_Upper -- ---------------------- function Sum_Strict_Upper (M : Square_Matrix) return Real is Sum : Real := 0.0; begin for Row in 1 .. N - 1 loop for Col in Row + 1 .. N loop Sum := Sum + abs M (Row, Col); end loop; end loop; return Sum; end Sum_Strict_Upper; M : Square_Matrix := A; -- Work space for solving eigensystem Threshold : Real; Sum : Real; Diag : Real_Vector (1 .. N); Diag_Adj : Real_Vector (1 .. N); -- The vector Diag_Adj indicates the amount of change in each value, -- while Diag tracks the value itself and Values holds the values as -- they were at the beginning. As the changes typically will be small -- compared to the absolute value of Diag, at the end of each iteration -- Diag is computed as Diag + Diag_Adj thus avoiding accumulating -- rounding errors. This technique is due to Rutishauser. begin if Compute_Vectors and then (Vectors'Length (1) /= N or else Vectors'Length (2) /= N) then raise Constraint_Error with "incompatible matrix dimensions"; elsif Values'Length /= N then raise Constraint_Error with "incompatible vector length"; elsif not Is_Symmetric (M) then raise Constraint_Error with "matrix not symmetric"; end if; -- Note: Only the locally declared matrix M and vectors (Diag, Diag_Adj) -- have lower bound equal to 1. The Vectors matrix may have -- different bounds, so take care indexing elements. Assignment -- as a whole is fine as sliding is automatic in that case. Vectors := (if not Compute_Vectors then (1 .. 0 => (1 .. 0 => 0.0)) else Unit_Matrix (Vectors'Length (1), Vectors'Length (2))); Values := Diagonal (M); Sweep : for Iteration in 1 .. Max_Iterations loop -- The first three iterations, perform rotation for any non-zero -- element. After this, rotate only for those that are not much -- smaller than the average off-diagnal element. After the fifth -- iteration, additionally zero out off-diagonal elements that are -- very small compared to elements on the diagonal with the same -- column or row index. Sum := Sum_Strict_Upper (M); exit Sweep when Sum = 0.0; Threshold := (if Iteration < 4 then 0.2 * Sum / Real (N**2) else 0.0); -- Iterate over all off-diagonal elements, rotating any that have -- an absolute value that exceeds the threshold. Diag := Values; Diag_Adj := (others => 0.0); -- Accumulates adjustments to Diag for Row in 1 .. N - 1 loop for Col in Row + 1 .. N loop -- If, before the rotation M (Row, Col) is tiny compared to -- Diag (Row) and Diag (Col), rotation is skipped. This is -- meaningful, as it produces no larger error than would be -- produced anyhow if the rotation had been performed. -- Suppress this optimization in the first four sweeps, so -- that this procedure can be used for computing eigenvectors -- of perturbed diagonal matrices. if Iteration > 4 and then Is_Tiny (M (Row, Col), Compared_To => Diag (Row)) and then Is_Tiny (M (Row, Col), Compared_To => Diag (Col)) then M (Row, Col) := 0.0; elsif abs M (Row, Col) > Threshold then Perform_Rotation : declare Tan : constant Real := Compute_Tan (M (Row, Col), Diag (Col) - Diag (Row)); Cos : constant Real := 1.0 / Sqrt (1.0 + Tan**2); Sin : constant Real := Tan * Cos; Tau : constant Real := Sin / (1.0 + Cos); Adj : constant Real := Tan * M (Row, Col); begin Diag_Adj (Row) := Diag_Adj (Row) - Adj; Diag_Adj (Col) := Diag_Adj (Col) + Adj; Diag (Row) := Diag (Row) - Adj; Diag (Col) := Diag (Col) + Adj; M (Row, Col) := 0.0; for J in 1 .. Row - 1 loop -- 1 <= J < Row Rotate (M (J, Row), M (J, Col), Sin, Tau); end loop; for J in Row + 1 .. Col - 1 loop -- Row < J < Col Rotate (M (Row, J), M (J, Col), Sin, Tau); end loop; for J in Col + 1 .. N loop -- Col < J <= N Rotate (M (Row, J), M (Col, J), Sin, Tau); end loop; for J in Vectors'Range (1) loop Rotate (Vectors (J, Row - 1 + Vectors'First (2)), Vectors (J, Col - 1 + Vectors'First (2)), Sin, Tau); end loop; end Perform_Rotation; end if; end loop; end loop; Values := Values + Diag_Adj; end loop Sweep; -- All normal matrices with valid values should converge perfectly. if Sum /= 0.0 then raise Constraint_Error with "eigensystem solution does not converge"; end if; end Jacobi; ----------- -- Solve -- ----------- function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector renames Instantiations.Solve; function Solve (A, X : Real_Matrix) return Real_Matrix renames Instantiations.Solve; ---------------------- -- Sort_Eigensystem -- ---------------------- procedure Sort_Eigensystem (Values : in out Real_Vector; Vectors : in out Real_Matrix) is procedure Swap (Left, Right : Integer); -- Swap Values (Left) with Values (Right), and also swap the -- corresponding eigenvectors. Note that lowerbounds may differ. function Less (Left, Right : Integer) return Boolean is (Values (Left) > Values (Right)); -- Sort by decreasing eigenvalue, see RM G.3.1 (76). procedure Sort is new Generic_Anonymous_Array_Sort (Integer); -- Sorts eigenvalues and eigenvectors by decreasing value procedure Swap (Left, Right : Integer) is begin Swap (Values (Left), Values (Right)); Swap_Column (Vectors, Left - Values'First + Vectors'First (2), Right - Values'First + Vectors'First (2)); end Swap; begin Sort (Values'First, Values'Last); end Sort_Eigensystem; --------------- -- Transpose -- --------------- function Transpose (X : Real_Matrix) return Real_Matrix is begin return R : Real_Matrix (X'Range (2), X'Range (1)) do Transpose (X, R); end return; end Transpose; ----------------- -- Unit_Matrix -- ----------------- function Unit_Matrix (Order : Positive; First_1 : Integer := 1; First_2 : Integer := 1) return Real_Matrix renames Instantiations.Unit_Matrix; ----------------- -- Unit_Vector -- ----------------- function Unit_Vector (Index : Integer; Order : Positive; First : Integer := 1) return Real_Vector renames Instantiations.Unit_Vector; end Ada.Numerics.Generic_Real_Arrays;