Mercurial > hg > CbC > CbC_gcc

comparison gcc/ada/libgnat/a-ngrear.adb @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents
children 84e7813d76e9
comparison
equal deleted inserted replaced
68:561a7518be6b 111:04ced10e8804
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-2017, 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;