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