| ------------------------------------------------------------------------------ |
| -- -- |
| -- GNAT RUN-TIME COMPONENTS -- |
| -- -- |
| -- ADA.NUMERICS.GENERIC_REAL_ARRAYS -- |
| -- -- |
| -- B o d y -- |
| -- -- |
| -- Copyright (C) 2006-2021, 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 -- |
| -- <http://www.gnu.org/licenses/>. -- |
| -- -- |
| -- 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))); |
| pragma Annotate |
| (CodePeer, False_Positive, "divide by zero", "H, P /= 0"); |
| |
| 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; |
