gcc/ada/libgnat/a-ngrear.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                     ADA.NUMERICS.GENERIC_REAL_ARRAYS                     --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --          Copyright (C) 2006-2023, Free Software Foundation, Inc.         --
  10 --                                                                          --
  11 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12 -- terms of the  GNU General Public License as published  by the Free Soft- --
  13 -- ware  Foundation;  either version 3,  or (at your option) any later ver- --
  14 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  15 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  16 -- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
  17 --                                                                          --
  18 -- As a special exception under Section 7 of GPL version 3, you are granted --
  19 -- additional permissions described in the GCC Runtime Library Exception,   --
  20 -- version 3.1, as published by the Free Software Foundation.               --
  21 --                                                                          --
  22 -- You should have received a copy of the GNU General Public License and    --
  23 -- a copy of the GCC Runtime Library Exception along with this program;     --
  24 -- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
  25 -- <http://www.gnu.org/licenses/>.                                          --
  26 --                                                                          --
  27 -- GNAT was originally developed  by the GNAT team at  New York University. --
  28 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  29 --                                                                          --
  30 ------------------------------------------------------------------------------
  31
  32 --  This version of Generic_Real_Arrays avoids the use of BLAS and LAPACK. One
  33 --  reason for this is new Ada 2012 requirements that prohibit algorithms such
  34 --  as Strassen's algorithm, which may be used by some BLAS implementations. In
  35 --  addition, some platforms lacked suitable compilers to compile the reference
  36 --  BLAS/LAPACK implementation. Finally, on some platforms there are more
  37 --  floating point types than supported by BLAS/LAPACK.
  38
  39 --  Preconditions, postconditions, ghost code, loop invariants and assertions
  40 --  in this unit are meant for analysis only, not for run-time checking, as it
  41 --  would be too costly otherwise. This is enforced by setting the assertion
  42 --  policy to Ignore.
  43
  44 pragma Assertion_Policy (Pre            => Ignore,
  45                          Post           => Ignore,
  46                          Ghost          => Ignore,
  47                          Loop_Invariant => Ignore,
  48                          Assert         => Ignore);
  49
  50 with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers;
  51
  52 with System; use System;
  53 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
  54
  55 package body Ada.Numerics.Generic_Real_Arrays is
  56
  57    package Ops renames System.Generic_Array_Operations;
  58
  59    function Is_Non_Zero (X : Real'Base) return Boolean is (X /= 0.0);
  60
  61    procedure Back_Substitute is new Ops.Back_Substitute
  62      (Scalar        => Real'Base,
  63       Matrix        => Real_Matrix,
  64       Is_Non_Zero   => Is_Non_Zero);
  65
  66    function Diagonal is new Ops.Diagonal
  67      (Scalar       => Real'Base,
  68       Vector       => Real_Vector,
  69       Matrix       => Real_Matrix);
  70
  71    procedure Forward_Eliminate is new Ops.Forward_Eliminate
  72     (Scalar        => Real'Base,
  73      Real          => Real'Base,
  74      Matrix        => Real_Matrix,
  75      Zero          => 0.0,
  76      One           => 1.0);
  77
  78    procedure Swap_Column is new Ops.Swap_Column
  79     (Scalar        => Real'Base,
  80      Matrix        => Real_Matrix);
  81
  82    procedure Transpose is new  Ops.Transpose
  83        (Scalar        => Real'Base,
  84         Matrix        => Real_Matrix);
  85
  86    function Is_Symmetric (A : Real_Matrix) return Boolean is
  87      (Transpose (A) = A);
  88    --  Return True iff A is symmetric, see RM G.3.1 (90).
  89
  90    function Is_Tiny (Value, Compared_To : Real) return Boolean is
  91      (abs Compared_To + 100.0 * abs (Value) = abs Compared_To);
  92    --  Return True iff the Value is much smaller in magnitude than the least
  93    --  significant digit of Compared_To.
  94
  95    procedure Jacobi
  96      (A               : Real_Matrix;
  97       Values          : out Real_Vector;
  98       Vectors         : out Real_Matrix;
  99       Compute_Vectors : Boolean := True);
 100    --  Perform Jacobi's eigensystem algorithm on real symmetric matrix A
 101
 102    function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
 103    --  Helper function that raises a Constraint_Error is the argument is
 104    --  not a square matrix, and otherwise returns its length.
 105
 106    procedure Rotate (X, Y : in out Real; Sin, Tau : Real);
 107    --  Perform a Givens rotation
 108
 109    procedure Sort_Eigensystem
 110      (Values  : in out Real_Vector;
 111       Vectors : in out Real_Matrix);
 112    --  Sort Values and associated Vectors by decreasing absolute value
 113
 114    procedure Swap (Left, Right : in out Real);
 115    --  Exchange Left and Right
 116
 117    function Sqrt is new Ops.Sqrt (Real);
 118    --  Instant a generic square root implementation here, in order to avoid
 119    --  instantiating a complete copy of Generic_Elementary_Functions.
 120    --  Speed of the square root is not a big concern here.
 121
 122    ------------
 123    -- Rotate --
 124    ------------
 125
 126    procedure Rotate (X, Y : in out Real; Sin, Tau : Real) is
 127       Old_X : constant Real := X;
 128       Old_Y : constant Real := Y;
 129    begin
 130       X := Old_X - Sin * (Old_Y + Old_X * Tau);
 131       Y := Old_Y + Sin * (Old_X - Old_Y * Tau);
 132    end Rotate;
 133
 134    ----------
 135    -- Swap --
 136    ----------
 137
 138    procedure Swap (Left, Right : in out Real) is
 139       Temp : constant Real := Left;
 140    begin
 141       Left := Right;
 142       Right := Temp;
 143    end Swap;
 144
 145    --  Instantiating the following subprograms directly would lead to
 146    --  name clashes, so use a local package.
 147
 148    package Instantiations is
 149
 150       function "+" is new
 151         Vector_Elementwise_Operation
 152           (X_Scalar      => Real'Base,
 153            Result_Scalar => Real'Base,
 154            X_Vector      => Real_Vector,
 155            Result_Vector => Real_Vector,
 156            Operation     => "+");
 157
 158       function "+" is new
 159         Matrix_Elementwise_Operation
 160           (X_Scalar      => Real'Base,
 161            Result_Scalar => Real'Base,
 162            X_Matrix      => Real_Matrix,
 163            Result_Matrix => Real_Matrix,
 164            Operation     => "+");
 165
 166       function "+" is new
 167         Vector_Vector_Elementwise_Operation
 168           (Left_Scalar   => Real'Base,
 169            Right_Scalar  => Real'Base,
 170            Result_Scalar => Real'Base,
 171            Left_Vector   => Real_Vector,
 172            Right_Vector  => Real_Vector,
 173            Result_Vector => Real_Vector,
 174            Operation     => "+");
 175
 176       function "+" is new
 177         Matrix_Matrix_Elementwise_Operation
 178           (Left_Scalar   => Real'Base,
 179            Right_Scalar  => Real'Base,
 180            Result_Scalar => Real'Base,
 181            Left_Matrix   => Real_Matrix,
 182            Right_Matrix  => Real_Matrix,
 183            Result_Matrix => Real_Matrix,
 184            Operation     => "+");
 185
 186       function "-" is new
 187         Vector_Elementwise_Operation
 188           (X_Scalar      => Real'Base,
 189            Result_Scalar => Real'Base,
 190            X_Vector      => Real_Vector,
 191            Result_Vector => Real_Vector,
 192            Operation     => "-");
 193
 194       function "-" is new
 195         Matrix_Elementwise_Operation
 196           (X_Scalar      => Real'Base,
 197            Result_Scalar => Real'Base,
 198            X_Matrix      => Real_Matrix,
 199            Result_Matrix => Real_Matrix,
 200            Operation     => "-");
 201
 202       function "-" is new
 203         Vector_Vector_Elementwise_Operation
 204           (Left_Scalar   => Real'Base,
 205            Right_Scalar  => Real'Base,
 206            Result_Scalar => Real'Base,
 207            Left_Vector   => Real_Vector,
 208            Right_Vector  => Real_Vector,
 209            Result_Vector => Real_Vector,
 210            Operation     => "-");
 211
 212       function "-" is new
 213         Matrix_Matrix_Elementwise_Operation
 214           (Left_Scalar   => Real'Base,
 215            Right_Scalar  => Real'Base,
 216            Result_Scalar => Real'Base,
 217            Left_Matrix   => Real_Matrix,
 218            Right_Matrix  => Real_Matrix,
 219            Result_Matrix => Real_Matrix,
 220            Operation     => "-");
 221
 222       function "*" is new
 223         Scalar_Vector_Elementwise_Operation
 224           (Left_Scalar   => Real'Base,
 225            Right_Scalar  => Real'Base,
 226            Result_Scalar => Real'Base,
 227            Right_Vector  => Real_Vector,
 228            Result_Vector => Real_Vector,
 229            Operation     => "*");
 230
 231       function "*" is new
 232         Scalar_Matrix_Elementwise_Operation
 233           (Left_Scalar   => Real'Base,
 234            Right_Scalar  => Real'Base,
 235            Result_Scalar => Real'Base,
 236            Right_Matrix  => Real_Matrix,
 237            Result_Matrix => Real_Matrix,
 238            Operation     => "*");
 239
 240       function "*" is new
 241         Vector_Scalar_Elementwise_Operation
 242           (Left_Scalar   => Real'Base,
 243            Right_Scalar  => Real'Base,
 244            Result_Scalar => Real'Base,
 245            Left_Vector   => Real_Vector,
 246            Result_Vector => Real_Vector,
 247            Operation     => "*");
 248
 249       function "*" is new
 250         Matrix_Scalar_Elementwise_Operation
 251           (Left_Scalar   => Real'Base,
 252            Right_Scalar  => Real'Base,
 253            Result_Scalar => Real'Base,
 254            Left_Matrix   => Real_Matrix,
 255            Result_Matrix => Real_Matrix,
 256            Operation     => "*");
 257
 258       function "*" is new
 259         Outer_Product
 260           (Left_Scalar   => Real'Base,
 261            Right_Scalar  => Real'Base,
 262            Result_Scalar => Real'Base,
 263            Left_Vector   => Real_Vector,
 264            Right_Vector  => Real_Vector,
 265            Matrix        => Real_Matrix);
 266
 267       function "*" is new
 268         Inner_Product
 269           (Left_Scalar   => Real'Base,
 270            Right_Scalar  => Real'Base,
 271            Result_Scalar => Real'Base,
 272            Left_Vector   => Real_Vector,
 273            Right_Vector  => Real_Vector,
 274            Zero          => 0.0);
 275
 276       function "*" is new
 277         Matrix_Vector_Product
 278           (Left_Scalar   => Real'Base,
 279            Right_Scalar  => Real'Base,
 280            Result_Scalar => Real'Base,
 281            Matrix        => Real_Matrix,
 282            Right_Vector  => Real_Vector,
 283            Result_Vector => Real_Vector,
 284            Zero          => 0.0);
 285
 286       function "*" is new
 287         Vector_Matrix_Product
 288           (Left_Scalar   => Real'Base,
 289            Right_Scalar  => Real'Base,
 290            Result_Scalar => Real'Base,
 291            Left_Vector   => Real_Vector,
 292            Matrix        => Real_Matrix,
 293            Result_Vector => Real_Vector,
 294            Zero          => 0.0);
 295
 296       function "*" is new
 297         Matrix_Matrix_Product
 298           (Left_Scalar   => Real'Base,
 299            Right_Scalar  => Real'Base,
 300            Result_Scalar => Real'Base,
 301            Left_Matrix   => Real_Matrix,
 302            Right_Matrix  => Real_Matrix,
 303            Result_Matrix => Real_Matrix,
 304            Zero          => 0.0);
 305
 306       function "/" is new
 307         Vector_Scalar_Elementwise_Operation
 308           (Left_Scalar   => Real'Base,
 309            Right_Scalar  => Real'Base,
 310            Result_Scalar => Real'Base,
 311            Left_Vector   => Real_Vector,
 312            Result_Vector => Real_Vector,
 313            Operation     => "/");
 314
 315       function "/" is new
 316         Matrix_Scalar_Elementwise_Operation
 317           (Left_Scalar   => Real'Base,
 318            Right_Scalar  => Real'Base,
 319            Result_Scalar => Real'Base,
 320            Left_Matrix   => Real_Matrix,
 321            Result_Matrix => Real_Matrix,
 322            Operation     => "/");
 323
 324       function "abs" is new
 325         L2_Norm
 326           (X_Scalar      => Real'Base,
 327            Result_Real   => Real'Base,
 328            X_Vector      => Real_Vector,
 329            "abs"         => "+");
 330       --  While the L2_Norm by definition uses the absolute values of the
 331       --  elements of X_Vector, for real values the subsequent squaring
 332       --  makes this unnecessary, so we substitute the "+" identity function
 333       --  instead.
 334
 335       function "abs" is new
 336         Vector_Elementwise_Operation
 337           (X_Scalar      => Real'Base,
 338            Result_Scalar => Real'Base,
 339            X_Vector      => Real_Vector,
 340            Result_Vector => Real_Vector,
 341            Operation     => "abs");
 342
 343       function "abs" is new
 344         Matrix_Elementwise_Operation
 345           (X_Scalar      => Real'Base,
 346            Result_Scalar => Real'Base,
 347            X_Matrix      => Real_Matrix,
 348            Result_Matrix => Real_Matrix,
 349            Operation     => "abs");
 350
 351       function Solve is new
 352         Matrix_Vector_Solution (Real'Base, 0.0, Real_Vector, Real_Matrix);
 353
 354       function Solve is new
 355         Matrix_Matrix_Solution (Real'Base, 0.0, Real_Matrix);
 356
 357       function Unit_Matrix is new
 358         Generic_Array_Operations.Unit_Matrix
 359           (Scalar        => Real'Base,
 360            Matrix        => Real_Matrix,
 361            Zero          => 0.0,
 362            One           => 1.0);
 363
 364       function Unit_Vector is new
 365         Generic_Array_Operations.Unit_Vector
 366           (Scalar        => Real'Base,
 367            Vector        => Real_Vector,
 368            Zero          => 0.0,
 369            One           => 1.0);
 370
 371    end Instantiations;
 372
 373    ---------
 374    -- "+" --
 375    ---------
 376
 377    function "+" (Right : Real_Vector) return Real_Vector
 378      renames Instantiations."+";
 379
 380    function "+" (Right : Real_Matrix) return Real_Matrix
 381      renames Instantiations."+";
 382
 383    function "+" (Left, Right : Real_Vector) return Real_Vector
 384      renames Instantiations."+";
 385
 386    function "+" (Left, Right : Real_Matrix) return Real_Matrix
 387      renames Instantiations."+";
 388
 389    ---------
 390    -- "-" --
 391    ---------
 392
 393    function "-" (Right : Real_Vector) return Real_Vector
 394      renames Instantiations."-";
 395
 396    function "-" (Right : Real_Matrix) return Real_Matrix
 397      renames Instantiations."-";
 398
 399    function "-" (Left, Right : Real_Vector) return Real_Vector
 400      renames Instantiations."-";
 401
 402    function "-" (Left, Right : Real_Matrix) return Real_Matrix
 403       renames Instantiations."-";
 404
 405    ---------
 406    -- "*" --
 407    ---------
 408
 409    --  Scalar multiplication
 410
 411    function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector
 412      renames Instantiations."*";
 413
 414    function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector
 415      renames Instantiations."*";
 416
 417    function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix
 418      renames Instantiations."*";
 419
 420    function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
 421      renames Instantiations."*";
 422
 423    --  Vector multiplication
 424
 425    function "*" (Left, Right : Real_Vector) return Real'Base
 426      renames Instantiations."*";
 427
 428    function "*" (Left, Right : Real_Vector) return Real_Matrix
 429      renames Instantiations."*";
 430
 431    function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector
 432      renames Instantiations."*";
 433
 434    function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector
 435      renames Instantiations."*";
 436
 437    --  Matrix Multiplication
 438
 439    function "*" (Left, Right : Real_Matrix) return Real_Matrix
 440      renames Instantiations."*";
 441
 442    ---------
 443    -- "/" --
 444    ---------
 445
 446    function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector
 447      renames Instantiations."/";
 448
 449    function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
 450      renames Instantiations."/";
 451
 452    -----------
 453    -- "abs" --
 454    -----------
 455
 456    function "abs" (Right : Real_Vector) return Real'Base
 457      renames Instantiations."abs";
 458
 459    function "abs" (Right : Real_Vector) return Real_Vector
 460      renames Instantiations."abs";
 461
 462    function "abs" (Right : Real_Matrix) return Real_Matrix
 463      renames Instantiations."abs";
 464
 465    -----------------
 466    -- Determinant --
 467    -----------------
 468
 469    function Determinant (A : Real_Matrix) return Real'Base is
 470       M : Real_Matrix := A;
 471       B : Real_Matrix (A'Range (1), 1 .. 0);
 472       R : Real'Base;
 473    begin
 474       Forward_Eliminate (M, B, R);
 475       return R;
 476    end Determinant;
 477
 478    -----------------
 479    -- Eigensystem --
 480    -----------------
 481
 482    procedure Eigensystem
 483      (A       : Real_Matrix;
 484       Values  : out Real_Vector;
 485       Vectors : out Real_Matrix)
 486    is
 487    begin
 488       Jacobi (A, Values, Vectors, Compute_Vectors => True);
 489       Sort_Eigensystem (Values, Vectors);
 490    end Eigensystem;
 491
 492    -----------------
 493    -- Eigenvalues --
 494    -----------------
 495
 496    function Eigenvalues (A : Real_Matrix) return Real_Vector is
 497    begin
 498       return Values : Real_Vector (A'Range (1)) do
 499          declare
 500             Vectors : Real_Matrix (1 .. 0, 1 .. 0);
 501          begin
 502             Jacobi (A, Values, Vectors, Compute_Vectors => False);
 503             Sort_Eigensystem (Values, Vectors);
 504          end;
 505       end return;
 506    end Eigenvalues;
 507
 508    -------------
 509    -- Inverse --
 510    -------------
 511
 512    function Inverse (A : Real_Matrix) return Real_Matrix is
 513      (Solve (A, Unit_Matrix (Length (A),
 514                              First_1 => A'First (2),
 515                              First_2 => A'First (1))));
 516
 517    ------------
 518    -- Jacobi --
 519    ------------
 520
 521    procedure Jacobi
 522      (A               : Real_Matrix;
 523       Values          : out Real_Vector;
 524       Vectors         : out Real_Matrix;
 525       Compute_Vectors : Boolean := True)
 526    is
 527       --  This subprogram uses Carl Gustav Jacob Jacobi's iterative method
 528       --  for computing eigenvalues and eigenvectors and is based on
 529       --  Rutishauser's implementation.
 530
 531       --  The given real symmetric matrix is transformed iteratively to
 532       --  diagonal form through a sequence of appropriately chosen elementary
 533       --  orthogonal transformations, called Jacobi rotations here.
 534
 535       --  The Jacobi method produces a systematic decrease of the sum of the
 536       --  squares of off-diagonal elements. Convergence to zero is quadratic,
 537       --  both for this implementation, as for the classic method that doesn't
 538       --  use row-wise scanning for pivot selection.
 539
 540       --  The numerical stability and accuracy of Jacobi's method make it the
 541       --  best choice here, even though for large matrices other methods will
 542       --  be significantly more efficient in both time and space.
 543
 544       --  While the eigensystem computations are absolutely foolproof for all
 545       --  real symmetric matrices, in presence of invalid values, or similar
 546       --  exceptional situations it might not. In such cases the results cannot
 547       --  be trusted and Constraint_Error is raised.
 548
 549       --  Note: this implementation needs temporary storage for 2 * N + N**2
 550       --        values of type Real.
 551
 552       Max_Iterations  : constant := 50;
 553       N               : constant Natural := Length (A);
 554
 555       subtype Square_Matrix is Real_Matrix (1 .. N, 1 .. N);
 556
 557       --  In order to annihilate the M (Row, Col) element, the
 558       --  rotation parameters Cos and Sin are computed as
 559       --  follows:
 560
 561       --    Theta = Cot (2.0 * Phi)
 562       --          = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col))
 563
 564       --  Then Tan (Phi) as the smaller root (in modulus) of
 565
 566       --    T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large)
 567
 568       function Compute_Tan (Theta : Real) return Real is
 569          (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta));
 570
 571       function Compute_Tan (P, H : Real) return Real is
 572          (if Is_Tiny (P, Compared_To => H) then P / H
 573           else Compute_Tan (Theta => H / (2.0 * P)));
 574       pragma Annotate
 575         (CodePeer, False_Positive, "divide by zero", "H, P /= 0");
 576
 577       function Sum_Strict_Upper (M : Square_Matrix) return Real;
 578       --  Return the sum of all elements in the strict upper triangle of M
 579
 580       ----------------------
 581       -- Sum_Strict_Upper --
 582       ----------------------
 583
 584       function Sum_Strict_Upper (M : Square_Matrix) return Real is
 585          Sum : Real := 0.0;
 586
 587       begin
 588          for Row in 1 .. N - 1 loop
 589             for Col in Row + 1 .. N loop
 590                Sum := Sum + abs M (Row, Col);
 591             end loop;
 592          end loop;
 593
 594          return Sum;
 595       end Sum_Strict_Upper;
 596
 597       M         : Square_Matrix := A; --  Work space for solving eigensystem
 598       Threshold : Real;
 599       Sum       : Real;
 600       Diag      : Real_Vector (1 .. N);
 601       Diag_Adj  : Real_Vector (1 .. N);
 602
 603       --  The vector Diag_Adj indicates the amount of change in each value,
 604       --  while Diag tracks the value itself and Values holds the values as
 605       --  they were at the beginning. As the changes typically will be small
 606       --  compared to the absolute value of Diag, at the end of each iteration
 607       --  Diag is computed as Diag + Diag_Adj thus avoiding accumulating
 608       --  rounding errors. This technique is due to Rutishauser.
 609
 610    begin
 611       if Compute_Vectors
 612          and then (Vectors'Length (1) /= N or else Vectors'Length (2) /= N)
 613       then
 614          raise Constraint_Error with "incompatible matrix dimensions";
 615
 616       elsif Values'Length /= N then
 617          raise Constraint_Error with "incompatible vector length";
 618
 619       elsif not Is_Symmetric (M) then
 620          raise Constraint_Error with "matrix not symmetric";
 621       end if;
 622
 623       --  Note: Only the locally declared matrix M and vectors (Diag, Diag_Adj)
 624       --        have lower bound equal to 1. The Vectors matrix may have
 625       --        different bounds, so take care indexing elements. Assignment
 626       --        as a whole is fine as sliding is automatic in that case.
 627
 628       Vectors := (if not Compute_Vectors then [1 .. 0 => [1 .. 0 => 0.0]]
 629                   else Unit_Matrix (Vectors'Length (1), Vectors'Length (2)));
 630       Values := Diagonal (M);
 631
 632       Sweep : for Iteration in 1 .. Max_Iterations loop
 633
 634          --  The first three iterations, perform rotation for any non-zero
 635          --  element. After this, rotate only for those that are not much
 636          --  smaller than the average off-diagnal element. After the fifth
 637          --  iteration, additionally zero out off-diagonal elements that are
 638          --  very small compared to elements on the diagonal with the same
 639          --  column or row index.
 640
 641          Sum := Sum_Strict_Upper (M);
 642
 643          exit Sweep when Sum = 0.0;
 644
 645          Threshold := (if Iteration < 4 then 0.2 * Sum / Real (N**2) else 0.0);
 646
 647          --  Iterate over all off-diagonal elements, rotating any that have
 648          --  an absolute value that exceeds the threshold.
 649
 650          Diag := Values;
 651          Diag_Adj := [others => 0.0]; -- Accumulates adjustments to Diag
 652
 653          for Row in 1 .. N - 1 loop
 654             for Col in Row + 1 .. N loop
 655
 656                --  If, before the rotation M (Row, Col) is tiny compared to
 657                --  Diag (Row) and Diag (Col), rotation is skipped. This is
 658                --  meaningful, as it produces no larger error than would be
 659                --  produced anyhow if the rotation had been performed.
 660                --  Suppress this optimization in the first four sweeps, so
 661                --  that this procedure can be used for computing eigenvectors
 662                --  of perturbed diagonal matrices.
 663
 664                if Iteration > 4
 665                   and then Is_Tiny (M (Row, Col), Compared_To => Diag (Row))
 666                   and then Is_Tiny (M (Row, Col), Compared_To => Diag (Col))
 667                then
 668                   M (Row, Col) := 0.0;
 669
 670                elsif abs M (Row, Col) > Threshold then
 671                   Perform_Rotation : declare
 672                      Tan : constant Real := Compute_Tan (M (Row, Col),
 673                                                Diag (Col) - Diag (Row));
 674                      Cos : constant Real := 1.0 / Sqrt (1.0 + Tan**2);
 675                      Sin : constant Real := Tan * Cos;
 676                      Tau : constant Real := Sin / (1.0 + Cos);
 677                      Adj : constant Real := Tan * M (Row, Col);
 678
 679                   begin
 680                      Diag_Adj (Row) := Diag_Adj (Row) - Adj;
 681                      Diag_Adj (Col) := Diag_Adj (Col) + Adj;
 682                      Diag (Row) := Diag (Row) - Adj;
 683                      Diag (Col) := Diag (Col) + Adj;
 684
 685                      M (Row, Col) := 0.0;
 686
 687                      for J in 1 .. Row - 1 loop        --  1 <= J < Row
 688                         Rotate (M (J, Row), M (J, Col), Sin, Tau);
 689                      end loop;
 690
 691                      for J in Row + 1 .. Col - 1 loop  --  Row < J < Col
 692                         Rotate (M (Row, J), M (J, Col), Sin, Tau);
 693                      end loop;
 694
 695                      for J in Col + 1 .. N loop        --  Col < J <= N
 696                         Rotate (M (Row, J), M (Col, J), Sin, Tau);
 697                      end loop;
 698
 699                      for J in Vectors'Range (1) loop
 700                         Rotate (Vectors (J, Row - 1 + Vectors'First (2)),
 701                                 Vectors (J, Col - 1 + Vectors'First (2)),
 702                                 Sin, Tau);
 703                      end loop;
 704                   end Perform_Rotation;
 705                end if;
 706             end loop;
 707          end loop;
 708
 709          Values := Values + Diag_Adj;
 710       end loop Sweep;
 711
 712       --  All normal matrices with valid values should converge perfectly.
 713
 714       if Sum /= 0.0 then
 715          raise Constraint_Error with "eigensystem solution does not converge";
 716       end if;
 717    end Jacobi;
 718
 719    -----------
 720    -- Solve --
 721    -----------
 722
 723    function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector
 724       renames Instantiations.Solve;
 725
 726    function Solve (A, X : Real_Matrix) return Real_Matrix
 727       renames Instantiations.Solve;
 728
 729    ----------------------
 730    -- Sort_Eigensystem --
 731    ----------------------
 732
 733    procedure Sort_Eigensystem
 734      (Values  : in out Real_Vector;
 735       Vectors : in out Real_Matrix)
 736    is
 737       procedure Swap (Left, Right : Integer);
 738       --  Swap Values (Left) with Values (Right), and also swap the
 739       --  corresponding eigenvectors. Note that lowerbounds may differ.
 740
 741       function Less (Left, Right : Integer) return Boolean is
 742         (Values (Left) > Values (Right));
 743       --  Sort by decreasing eigenvalue, see RM G.3.1 (76).
 744
 745       procedure Sort is new Generic_Anonymous_Array_Sort (Integer);
 746       --  Sorts eigenvalues and eigenvectors by decreasing value
 747
 748       procedure Swap (Left, Right : Integer) is
 749       begin
 750          Swap (Values (Left), Values (Right));
 751          Swap_Column (Vectors, Left - Values'First + Vectors'First (2),
 752                                Right - Values'First + Vectors'First (2));
 753       end Swap;
 754
 755    begin
 756       Sort (Values'First, Values'Last);
 757    end Sort_Eigensystem;
 758
 759    ---------------
 760    -- Transpose --
 761    ---------------
 762
 763    function Transpose (X : Real_Matrix) return Real_Matrix is
 764    begin
 765       return R : Real_Matrix (X'Range (2), X'Range (1)) do
 766          Transpose (X, R);
 767       end return;
 768    end Transpose;
 769
 770    -----------------
 771    -- Unit_Matrix --
 772    -----------------
 773
 774    function Unit_Matrix
 775      (Order   : Positive;
 776       First_1 : Integer := 1;
 777       First_2 : Integer := 1) return Real_Matrix
 778      renames Instantiations.Unit_Matrix;
 779
 780    -----------------
 781    -- Unit_Vector --
 782    -----------------
 783
 784    function Unit_Vector
 785      (Index : Integer;
 786       Order : Positive;
 787       First : Integer := 1) return Real_Vector
 788      renames Instantiations.Unit_Vector;
 789
 790 end Ada.Numerics.Generic_Real_Arrays;