775 lines
27 KiB
Ada
775 lines
27 KiB
Ada
------------------------------------------------------------------------------
|
|
-- --
|
|
-- GNAT RUN-TIME COMPONENTS --
|
|
-- --
|
|
-- ADA.NUMERICS.GENERIC_REAL_ARRAYS --
|
|
-- --
|
|
-- B o d y --
|
|
-- --
|
|
-- Copyright (C) 2006-2012, 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, Real_Vector, Real_Matrix);
|
|
|
|
function Solve is new Matrix_Matrix_Solution (Real'Base, 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))));
|
|
|
|
------------
|
|
-- 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;
|