https://gcc.gnu.org/g:b65f4fe10996c247751b0636773009659d92ab6c
commit r12-11016-gb65f4fe10996c247751b0636773009659d92ab6c
Author: Eric Botcazou <ebotca...@adacore.com>
Date: Fri Apr 4 11:45:23 2025 +0200
Ada: Fix thinko in Eigensystem for complex Hermitian matrices
The implementation solves the eigensystem for a NxN complex Hermitian matrix
by first solving it for a 2Nx2N real symmetric matrix and then interpreting
the 2Nx1 real vectors as Nx1 complex ones, but the last step does not work.
The patch fixes the last step and also performs a small cleanup throughout
the implementation, mostly in the commentary and without functional changes.
gcc/ada/
* libgnat/a-ngcoar.adb (Eigensystem): Adjust notation and fix the
layout of the real symmetric matrix in the main comment. Adjust
the layout of the associated code accordingly and correctly turn
the 2Nx1 real vectors into Nx1 complex ones.
(Eigenvalues): Minor similar tweaks.
* libgnat/a-ngrear.adb (Jacobi): Minor tweaks in the main comment.
Adjust notation and corresponding parameter names of functions.
Fix call to Unit_Matrix routine. Adjust the comment describing
the various kinds of iterations to match the implementation.
Diff:
---
gcc/ada/libgnat/a-ngcoar.adb | 42 +++++++++++-----------
gcc/ada/libgnat/a-ngrear.adb | 85 ++++++++++++++++++++++++--------------------
2 files changed, 66 insertions(+), 61 deletions(-)
diff --git a/gcc/ada/libgnat/a-ngcoar.adb b/gcc/ada/libgnat/a-ngcoar.adb
index 8dfbc3b174a1..cbbc6d274aa3 100644
--- a/gcc/ada/libgnat/a-ngcoar.adb
+++ b/gcc/ada/libgnat/a-ngcoar.adb
@@ -1058,19 +1058,21 @@ package body Ada.Numerics.Generic_Complex_Arrays is
is
N : constant Natural := Length (A);
- -- For a Hermitian matrix C, we convert the eigenvalue problem to a
- -- real symmetric one: if C = A + i * B, then the (N, N) complex
+ -- For a Hermitian matrix A, we convert the eigenvalue problem to a
+ -- real symmetric one: if A = X + i * Y, then the (N, N) complex
-- eigenvalue problem:
- -- (A + i * B) * (u + i * v) = Lambda * (u + i * v)
+ --
+ -- (X + i * Y) * (u + i * v) = Lambda * (u + i * v)
--
-- is equivalent to the (2 * N, 2 * N) real eigenvalue problem:
- -- [ A, B ] [ u ] = Lambda * [ u ]
- -- [ -B, A ] [ v ] [ v ]
--
- -- Note that the (2 * N, 2 * N) matrix above is symmetric, as
- -- Transpose (A) = A and Transpose (B) = -B if C is Hermitian.
+ -- [ X, -Y ] [ u ] = Lambda * [ u ]
+ -- [ Y, X ] [ v ] [ v ]
+ --
+ -- Note that the (2 * N, 2 * N) matrix M above is symmetric, because
+ -- Transpose (X) = X and Transpose (Y) = -Y as A is Hermitian.
- -- We solve this eigensystem using the real-valued algorithms. The final
+ -- We solve this eigensystem using the real-valued algorithm. The final
-- result will have every eigenvalue twice, so in the sorted output we
-- just pick every second value, with associated eigenvector u + i * v.
@@ -1085,10 +1087,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
C : constant Complex :=
(A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
begin
- M (J, K) := Re (C);
- M (J + N, K + N) := Re (C);
- M (J + N, K) := Im (C);
- M (J, K + N) := -Im (C);
+ M (J, K) := Re (C); M (J, K + N) := -Im (C);
+ M (J + N, K) := Im (C); M (J + N, K + N) := Re (C);
end;
end loop;
end loop;
@@ -1103,10 +1103,9 @@ package body Ada.Numerics.Generic_Complex_Arrays is
for K in 1 .. N loop
declare
- Row : constant Integer := Vectors'First (2) + (K - 1);
+ Row : constant Integer := Vectors'First (1) + (K - 1);
begin
- Vectors (Row, Col)
- := (Vecs (J * 2, Col), Vecs (J * 2, Col + N));
+ Vectors (Row, Col) := (Vecs (K, 2 * J), Vecs (K + N, 2 * J));
end;
end loop;
end;
@@ -1118,13 +1117,14 @@ package body Ada.Numerics.Generic_Complex_Arrays is
-----------------
function Eigenvalues (A : Complex_Matrix) return Real_Vector is
- -- See Eigensystem for a description of the algorithm
-
N : constant Natural := Length (A);
- R : Real_Vector (A'Range (1));
+
+ -- See Eigensystem for a description of the algorithm
M : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
+ R : Real_Vector (A'Range (1));
Vals : Real_Vector (1 .. 2 * N);
+
begin
for J in 1 .. N loop
for K in 1 .. N loop
@@ -1132,10 +1132,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
C : constant Complex :=
(A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
begin
- M (J, K) := Re (C);
- M (J + N, K + N) := Re (C);
- M (J + N, K) := Im (C);
- M (J, K + N) := -Im (C);
+ M (J, K) := Re (C); M (J, K + N) := -Im (C);
+ M (J + N, K) := Im (C); M (J + N, K + N) := Re (C);
end;
end loop;
end loop;
diff --git a/gcc/ada/libgnat/a-ngrear.adb b/gcc/ada/libgnat/a-ngrear.adb
index 844d6264ee72..1d07a8a9b5fe 100644
--- a/gcc/ada/libgnat/a-ngrear.adb
+++ b/gcc/ada/libgnat/a-ngrear.adb
@@ -85,7 +85,7 @@ package body Ada.Numerics.Generic_Real_Arrays is
function Is_Symmetric (A : Real_Matrix) return Boolean is
(Transpose (A) = A);
- -- Return True iff A is symmetric, see RM G.3.1 (90).
+ -- 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);
@@ -104,7 +104,7 @@ package body Ada.Numerics.Generic_Real_Arrays is
-- not a square matrix, and otherwise returns its length.
procedure Rotate (X, Y : in out Real; Sin, Tau : Real);
- -- Perform a Givens rotation
+ -- Perform a Givens rotation of angle Theta given by sin and sin/(1 + cos)
procedure Sort_Eigensystem
(Values : in out Real_Vector;
@@ -525,13 +525,13 @@ package body Ada.Numerics.Generic_Real_Arrays is
Vectors : out Real_Matrix;
Compute_Vectors : Boolean)
is
- -- This subprogram uses Carl Gustav Jacob Jacobi's iterative method
- -- for computing eigenvalues and eigenvectors and is based on
+ -- This subprogram uses Carl Gustav Jacob Jacobi's cyclic 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.
+ -- orthogonal transformations, called Jacobi rotations.
-- The Jacobi method produces a systematic decrease of the sum of the
-- squares of off-diagonal elements. Convergence to zero is quadratic,
@@ -542,41 +542,44 @@ package body Ada.Numerics.Generic_Real_Arrays is
-- 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
+ -- While the eigensystem computation is 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.
+ -- exceptional situations, it may not be. 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);
+ 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:
+ -- In order to annihilate the M (Row, Col) element, the rotation angle
+ -- Theta is chosen as follows:
- -- Theta = Cot (2.0 * Phi)
- -- = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col))
+ -- Cot (2.0 * Theta) = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col))
- -- Then Tan (Phi) as the smaller root (in modulus) of
+ -- If C = Cot (2.0 * Theta), then Tan (Theta) is computed as the smaller
+ -- root (in modulus) of:
- -- T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large)
+ -- X**2 + 2 * C * X - 1 = 0
- function Compute_Tan (Theta : Real) return Real is
- (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta));
+ -- or else as 0.5 / C, if C is large.
+
+ function Compute_Tan (C : Real) return Real is
+ (Real'Copy_Sign (1.0 / (abs C + Sqrt (1.0 + C**2)), C));
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)));
+ (if Is_Tiny (P, Compared_To => H)
+ then P / H
+ else Compute_Tan (C => 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
+ -- Return the sum of the absolute value of all the elements in the
+ -- strict upper triangle of M.
----------------------
-- Sum_Strict_Upper --
@@ -595,11 +598,13 @@ package body Ada.Numerics.Generic_Real_Arrays is
return Sum;
end Sum_Strict_Upper;
+ -- Local variables
+
M : Square_Matrix := A; -- Work space for solving eigensystem
- Threshold : Real;
- Sum : Real;
Diag : Real_Vector (1 .. N);
Diag_Adj : Real_Vector (1 .. N);
+ Sum : Real;
+ Threshold : Real;
-- The vector Diag_Adj indicates the amount of change in each value,
-- while Diag tracks the value itself and Values holds the values as
@@ -621,22 +626,24 @@ package body Ada.Numerics.Generic_Real_Arrays is
raise Constraint_Error with "matrix not symmetric";
end if;
+ Values := Diagonal (M);
+
-- 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);
+ Vectors := (if Compute_Vectors
+ then Unit_Matrix (N)
+ else [1 .. 0 => [1 .. 0 => 0.0]]);
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
+ -- During the first three iterations, perform the rotation only for
+ -- elements that are not much smaller than the average off-diagonal
+ -- element. After this, rotate for any non-zero elements. 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);
@@ -645,8 +652,8 @@ package body Ada.Numerics.Generic_Real_Arrays is
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.
+ -- 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
@@ -654,11 +661,11 @@ package body Ada.Numerics.Generic_Real_Arrays is
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
+ -- 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
+ -- Suppress this optimization in the first four iterations, so
-- that this procedure can be used for computing eigenvectors
-- of perturbed diagonal matrices.
@@ -670,8 +677,8 @@ package body Ada.Numerics.Generic_Real_Arrays is
elsif abs M (Row, Col) > Threshold then
Perform_Rotation : declare
- Tan : constant Real := Compute_Tan (M (Row, Col),
- Diag (Col) - Diag (Row));
+ 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);
@@ -710,7 +717,7 @@ package body Ada.Numerics.Generic_Real_Arrays is
Values := Values + Diag_Adj;
end loop Sweep;
- -- All normal matrices with valid values should converge perfectly.
+ -- All normal matrices with valid values should converge perfectly
if Sum /= 0.0 then
raise Constraint_Error with "eigensystem solution does not converge";
