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-2012, 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 341 new Matrix_Vector_Solution (Real'Base, Real_Vector, Real_Matrix); 342 343 function Solve is new Matrix_Matrix_Solution (Real'Base, Real_Matrix); 344 345 function Unit_Matrix is new 346 Generic_Array_Operations.Unit_Matrix 347 (Scalar => Real'Base, 348 Matrix => Real_Matrix, 349 Zero => 0.0, 350 One => 1.0); 351 352 function Unit_Vector is new 353 Generic_Array_Operations.Unit_Vector 354 (Scalar => Real'Base, 355 Vector => Real_Vector, 356 Zero => 0.0, 357 One => 1.0); 358 359 end Instantiations; 360 361 --------- 362 -- "+" -- 363 --------- 364 365 function "+" (Right : Real_Vector) return Real_Vector 366 renames Instantiations."+"; 367 368 function "+" (Right : Real_Matrix) return Real_Matrix 369 renames Instantiations."+"; 370 371 function "+" (Left, Right : Real_Vector) return Real_Vector 372 renames Instantiations."+"; 373 374 function "+" (Left, Right : Real_Matrix) return Real_Matrix 375 renames Instantiations."+"; 376 377 --------- 378 -- "-" -- 379 --------- 380 381 function "-" (Right : Real_Vector) return Real_Vector 382 renames Instantiations."-"; 383 384 function "-" (Right : Real_Matrix) return Real_Matrix 385 renames Instantiations."-"; 386 387 function "-" (Left, Right : Real_Vector) return Real_Vector 388 renames Instantiations."-"; 389 390 function "-" (Left, Right : Real_Matrix) return Real_Matrix 391 renames Instantiations."-"; 392 393 --------- 394 -- "*" -- 395 --------- 396 397 -- Scalar multiplication 398 399 function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector 400 renames Instantiations."*"; 401 402 function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector 403 renames Instantiations."*"; 404 405 function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix 406 renames Instantiations."*"; 407 408 function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix 409 renames Instantiations."*"; 410 411 -- Vector multiplication 412 413 function "*" (Left, Right : Real_Vector) return Real'Base 414 renames Instantiations."*"; 415 416 function "*" (Left, Right : Real_Vector) return Real_Matrix 417 renames Instantiations."*"; 418 419 function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector 420 renames Instantiations."*"; 421 422 function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector 423 renames Instantiations."*"; 424 425 -- Matrix Multiplication 426 427 function "*" (Left, Right : Real_Matrix) return Real_Matrix 428 renames Instantiations."*"; 429 430 --------- 431 -- "/" -- 432 --------- 433 434 function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector 435 renames Instantiations."/"; 436 437 function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix 438 renames Instantiations."/"; 439 440 ----------- 441 -- "abs" -- 442 ----------- 443 444 function "abs" (Right : Real_Vector) return Real'Base 445 renames Instantiations."abs"; 446 447 function "abs" (Right : Real_Vector) return Real_Vector 448 renames Instantiations."abs"; 449 450 function "abs" (Right : Real_Matrix) return Real_Matrix 451 renames Instantiations."abs"; 452 453 ----------------- 454 -- Determinant -- 455 ----------------- 456 457 function Determinant (A : Real_Matrix) return Real'Base is 458 M : Real_Matrix := A; 459 B : Real_Matrix (A'Range (1), 1 .. 0); 460 R : Real'Base; 461 begin 462 Forward_Eliminate (M, B, R); 463 return R; 464 end Determinant; 465 466 ----------------- 467 -- Eigensystem -- 468 ----------------- 469 470 procedure Eigensystem 471 (A : Real_Matrix; 472 Values : out Real_Vector; 473 Vectors : out Real_Matrix) 474 is 475 begin 476 Jacobi (A, Values, Vectors, Compute_Vectors => True); 477 Sort_Eigensystem (Values, Vectors); 478 end Eigensystem; 479 480 ----------------- 481 -- Eigenvalues -- 482 ----------------- 483 484 function Eigenvalues (A : Real_Matrix) return Real_Vector is 485 begin 486 return Values : Real_Vector (A'Range (1)) do 487 declare 488 Vectors : Real_Matrix (1 .. 0, 1 .. 0); 489 begin 490 Jacobi (A, Values, Vectors, Compute_Vectors => False); 491 Sort_Eigensystem (Values, Vectors); 492 end; 493 end return; 494 end Eigenvalues; 495 496 ------------- 497 -- Inverse -- 498 ------------- 499 500 function Inverse (A : Real_Matrix) return Real_Matrix is 501 (Solve (A, Unit_Matrix (Length (A)))); 502 503 ------------ 504 -- Jacobi -- 505 ------------ 506 507 procedure Jacobi 508 (A : Real_Matrix; 509 Values : out Real_Vector; 510 Vectors : out Real_Matrix; 511 Compute_Vectors : Boolean := True) 512 is 513 -- This subprogram uses Carl Gustav Jacob Jacobi's iterative method 514 -- for computing eigenvalues and eigenvectors and is based on 515 -- Rutishauser's implementation. 516 517 -- The given real symmetric matrix is transformed iteratively to 518 -- diagonal form through a sequence of appropriately chosen elementary 519 -- orthogonal transformations, called Jacobi rotations here. 520 521 -- The Jacobi method produces a systematic decrease of the sum of the 522 -- squares of off-diagonal elements. Convergence to zero is quadratic, 523 -- both for this implementation, as for the classic method that doesn't 524 -- use row-wise scanning for pivot selection. 525 526 -- The numerical stability and accuracy of Jacobi's method make it the 527 -- best choice here, even though for large matrices other methods will 528 -- be significantly more efficient in both time and space. 529 530 -- While the eigensystem computations are absolutely foolproof for all 531 -- real symmetric matrices, in presence of invalid values, or similar 532 -- exceptional situations it might not. In such cases the results cannot 533 -- be trusted and Constraint_Error is raised. 534 535 -- Note: this implementation needs temporary storage for 2 * N + N**2 536 -- values of type Real. 537 538 Max_Iterations : constant := 50; 539 N : constant Natural := Length (A); 540 541 subtype Square_Matrix is Real_Matrix (1 .. N, 1 .. N); 542 543 -- In order to annihilate the M (Row, Col) element, the 544 -- rotation parameters Cos and Sin are computed as 545 -- follows: 546 547 -- Theta = Cot (2.0 * Phi) 548 -- = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col)) 549 550 -- Then Tan (Phi) as the smaller root (in modulus) of 551 552 -- T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large) 553 554 function Compute_Tan (Theta : Real) return Real is 555 (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta)); 556 557 function Compute_Tan (P, H : Real) return Real is 558 (if Is_Tiny (P, Compared_To => H) then P / H 559 else Compute_Tan (Theta => H / (2.0 * P))); 560 561 function Sum_Strict_Upper (M : Square_Matrix) return Real; 562 -- Return the sum of all elements in the strict upper triangle of M 563 564 ---------------------- 565 -- Sum_Strict_Upper -- 566 ---------------------- 567 568 function Sum_Strict_Upper (M : Square_Matrix) return Real is 569 Sum : Real := 0.0; 570 571 begin 572 for Row in 1 .. N - 1 loop 573 for Col in Row + 1 .. N loop 574 Sum := Sum + abs M (Row, Col); 575 end loop; 576 end loop; 577 578 return Sum; 579 end Sum_Strict_Upper; 580 581 M : Square_Matrix := A; -- Work space for solving eigensystem 582 Threshold : Real; 583 Sum : Real; 584 Diag : Real_Vector (1 .. N); 585 Diag_Adj : Real_Vector (1 .. N); 586 587 -- The vector Diag_Adj indicates the amount of change in each value, 588 -- while Diag tracks the value itself and Values holds the values as 589 -- they were at the beginning. As the changes typically will be small 590 -- compared to the absolute value of Diag, at the end of each iteration 591 -- Diag is computed as Diag + Diag_Adj thus avoiding accumulating 592 -- rounding errors. This technique is due to Rutishauser. 593 594 begin 595 if Compute_Vectors 596 and then (Vectors'Length (1) /= N or else Vectors'Length (2) /= N) 597 then 598 raise Constraint_Error with "incompatible matrix dimensions"; 599 600 elsif Values'Length /= N then 601 raise Constraint_Error with "incompatible vector length"; 602 603 elsif not Is_Symmetric (M) then 604 raise Constraint_Error with "matrix not symmetric"; 605 end if; 606 607 -- Note: Only the locally declared matrix M and vectors (Diag, Diag_Adj) 608 -- have lower bound equal to 1. The Vectors matrix may have 609 -- different bounds, so take care indexing elements. Assignment 610 -- as a whole is fine as sliding is automatic in that case. 611 612 Vectors := (if not Compute_Vectors then (1 .. 0 => (1 .. 0 => 0.0)) 613 else Unit_Matrix (Vectors'Length (1), Vectors'Length (2))); 614 Values := Diagonal (M); 615 616 Sweep : for Iteration in 1 .. Max_Iterations loop 617 618 -- The first three iterations, perform rotation for any non-zero 619 -- element. After this, rotate only for those that are not much 620 -- smaller than the average off-diagnal element. After the fifth 621 -- iteration, additionally zero out off-diagonal elements that are 622 -- very small compared to elements on the diagonal with the same 623 -- column or row index. 624 625 Sum := Sum_Strict_Upper (M); 626 627 exit Sweep when Sum = 0.0; 628 629 Threshold := (if Iteration < 4 then 0.2 * Sum / Real (N**2) else 0.0); 630 631 -- Iterate over all off-diagonal elements, rotating any that have 632 -- an absolute value that exceeds the threshold. 633 634 Diag := Values; 635 Diag_Adj := (others => 0.0); -- Accumulates adjustments to Diag 636 637 for Row in 1 .. N - 1 loop 638 for Col in Row + 1 .. N loop 639 640 -- If, before the rotation M (Row, Col) is tiny compared to 641 -- Diag (Row) and Diag (Col), rotation is skipped. This is 642 -- meaningful, as it produces no larger error than would be 643 -- produced anyhow if the rotation had been performed. 644 -- Suppress this optimization in the first four sweeps, so 645 -- that this procedure can be used for computing eigenvectors 646 -- of perturbed diagonal matrices. 647 648 if Iteration > 4 649 and then Is_Tiny (M (Row, Col), Compared_To => Diag (Row)) 650 and then Is_Tiny (M (Row, Col), Compared_To => Diag (Col)) 651 then 652 M (Row, Col) := 0.0; 653 654 elsif abs M (Row, Col) > Threshold then 655 Perform_Rotation : declare 656 Tan : constant Real := Compute_Tan (M (Row, Col), 657 Diag (Col) - Diag (Row)); 658 Cos : constant Real := 1.0 / Sqrt (1.0 + Tan**2); 659 Sin : constant Real := Tan * Cos; 660 Tau : constant Real := Sin / (1.0 + Cos); 661 Adj : constant Real := Tan * M (Row, Col); 662 663 begin 664 Diag_Adj (Row) := Diag_Adj (Row) - Adj; 665 Diag_Adj (Col) := Diag_Adj (Col) + Adj; 666 Diag (Row) := Diag (Row) - Adj; 667 Diag (Col) := Diag (Col) + Adj; 668 669 M (Row, Col) := 0.0; 670 671 for J in 1 .. Row - 1 loop -- 1 <= J < Row 672 Rotate (M (J, Row), M (J, Col), Sin, Tau); 673 end loop; 674 675 for J in Row + 1 .. Col - 1 loop -- Row < J < Col 676 Rotate (M (Row, J), M (J, Col), Sin, Tau); 677 end loop; 678 679 for J in Col + 1 .. N loop -- Col < J <= N 680 Rotate (M (Row, J), M (Col, J), Sin, Tau); 681 end loop; 682 683 for J in Vectors'Range (1) loop 684 Rotate (Vectors (J, Row - 1 + Vectors'First (2)), 685 Vectors (J, Col - 1 + Vectors'First (2)), 686 Sin, Tau); 687 end loop; 688 end Perform_Rotation; 689 end if; 690 end loop; 691 end loop; 692 693 Values := Values + Diag_Adj; 694 end loop Sweep; 695 696 -- All normal matrices with valid values should converge perfectly. 697 698 if Sum /= 0.0 then 699 raise Constraint_Error with "eigensystem solution does not converge"; 700 end if; 701 end Jacobi; 702 703 ----------- 704 -- Solve -- 705 ----------- 706 707 function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector 708 renames Instantiations.Solve; 709 710 function Solve (A, X : Real_Matrix) return Real_Matrix 711 renames Instantiations.Solve; 712 713 ---------------------- 714 -- Sort_Eigensystem -- 715 ---------------------- 716 717 procedure Sort_Eigensystem 718 (Values : in out Real_Vector; 719 Vectors : in out Real_Matrix) 720 is 721 procedure Swap (Left, Right : Integer); 722 -- Swap Values (Left) with Values (Right), and also swap the 723 -- corresponding eigenvectors. Note that lowerbounds may differ. 724 725 function Less (Left, Right : Integer) return Boolean is 726 (Values (Left) > Values (Right)); 727 -- Sort by decreasing eigenvalue, see RM G.3.1 (76). 728 729 procedure Sort is new Generic_Anonymous_Array_Sort (Integer); 730 -- Sorts eigenvalues and eigenvectors by decreasing value 731 732 procedure Swap (Left, Right : Integer) is 733 begin 734 Swap (Values (Left), Values (Right)); 735 Swap_Column (Vectors, Left - Values'First + Vectors'First (2), 736 Right - Values'First + Vectors'First (2)); 737 end Swap; 738 739 begin 740 Sort (Values'First, Values'Last); 741 end Sort_Eigensystem; 742 743 --------------- 744 -- Transpose -- 745 --------------- 746 747 function Transpose (X : Real_Matrix) return Real_Matrix is 748 begin 749 return R : Real_Matrix (X'Range (2), X'Range (1)) do 750 Transpose (X, R); 751 end return; 752 end Transpose; 753 754 ----------------- 755 -- Unit_Matrix -- 756 ----------------- 757 758 function Unit_Matrix 759 (Order : Positive; 760 First_1 : Integer := 1; 761 First_2 : Integer := 1) return Real_Matrix 762 renames Instantiations.Unit_Matrix; 763 764 ----------------- 765 -- Unit_Vector -- 766 ----------------- 767 768 function Unit_Vector 769 (Index : Integer; 770 Order : Positive; 771 First : Integer := 1) return Real_Vector 772 renames Instantiations.Unit_Vector; 773 774end Ada.Numerics.Generic_Real_Arrays; 775