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-2018, 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 39with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers; 40 41with System; use System; 42with System.Generic_Array_Operations; use System.Generic_Array_Operations; 43 44package 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 777end Ada.Numerics.Generic_Real_Arrays; 778