OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [a-ngrear.adb] - Blame information for rev 706

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 706 jeremybenn
------------------------------------------------------------------------------
2
--                                                                          --
3
--                         GNAT RUN-TIME COMPONENTS                         --
4
--                                                                          --
5
--                     ADA.NUMERICS.GENERIC_REAL_ARRAYS                     --
6
--                                                                          --
7
--                                 B o d y                                  --
8
--                                                                          --
9
--          Copyright (C) 2006-2011, 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
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
      Values  : Real_Vector (A'Range (1));
486
      Vectors : Real_Matrix (1 .. 0, 1 .. 0);
487
   begin
488
      Jacobi (A, Values, Vectors, Compute_Vectors => False);
489
      Sort_Eigensystem (Values, Vectors);
490
      return Values;
491
   end Eigenvalues;
492
 
493
   -------------
494
   -- Inverse --
495
   -------------
496
 
497
   function Inverse (A : Real_Matrix) return Real_Matrix is
498
     (Solve (A, Unit_Matrix (Length (A))));
499
 
500
   ------------
501
   -- Jacobi --
502
   ------------
503
 
504
   procedure Jacobi
505
     (A               : Real_Matrix;
506
      Values          : out Real_Vector;
507
      Vectors         : out Real_Matrix;
508
      Compute_Vectors : Boolean := True)
509
   is
510
      --  This subprogram uses Carl Gustav Jacob Jacobi's iterative method
511
      --  for computing eigenvalues and eigenvectors and is based on
512
      --  Rutishauser's implementation.
513
 
514
      --  The given real symmetric matrix is transformed iteratively to
515
      --  diagonal form through a sequence of appropriately chosen elementary
516
      --  orthogonal transformations, called Jacobi rotations here.
517
 
518
      --  The Jacobi method produces a systematic decrease of the sum of the
519
      --  squares of off-diagonal elements. Convergence to zero is quadratic,
520
      --  both for this implementation, as for the classic method that doesn't
521
      --  use row-wise scanning for pivot selection.
522
 
523
      --  The numerical stability and accuracy of Jacobi's method make it the
524
      --  best choice here, even though for large matrices other methods will
525
      --  be significantly more efficient in both time and space.
526
 
527
      --  While the eigensystem computations are absolutely foolproof for all
528
      --  real symmetric matrices, in presence of invalid values, or similar
529
      --  exceptional situations it might not. In such cases the results cannot
530
      --  be trusted and Constraint_Error is raised.
531
 
532
      --  Note: this implementation needs temporary storage for 2 * N + N**2
533
      --        values of type Real.
534
 
535
      Max_Iterations  : constant := 50;
536
      N               : constant Natural := Length (A);
537
 
538
      subtype Square_Matrix is Real_Matrix (1 .. N, 1 .. N);
539
 
540
      --  In order to annihilate the M (Row, Col) element, the
541
      --  rotation parameters Cos and Sin are computed as
542
      --  follows:
543
 
544
      --    Theta = Cot (2.0 * Phi)
545
      --          = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col))
546
 
547
      --  Then Tan (Phi) as the smaller root (in modulus) of
548
 
549
      --    T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large)
550
 
551
      function Compute_Tan (Theta : Real) return Real is
552
         (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta));
553
 
554
      function Compute_Tan (P, H : Real) return Real is
555
         (if Is_Tiny (P, Compared_To => H) then P / H
556
          else Compute_Tan (Theta => H / (2.0 * P)));
557
 
558
      function Sum_Strict_Upper (M : Square_Matrix) return Real;
559
      --  Return the sum of all elements in the strict upper triangle of M
560
 
561
      ----------------------
562
      -- Sum_Strict_Upper --
563
      ----------------------
564
 
565
      function Sum_Strict_Upper (M : Square_Matrix) return Real is
566
         Sum : Real := 0.0;
567
 
568
      begin
569
         for Row in 1 .. N - 1 loop
570
            for Col in Row + 1 .. N loop
571
               Sum := Sum + abs M (Row, Col);
572
            end loop;
573
         end loop;
574
 
575
         return Sum;
576
      end Sum_Strict_Upper;
577
 
578
      M         : Square_Matrix := A; --  Work space for solving eigensystem
579
      Threshold : Real;
580
      Sum       : Real;
581
      Diag      : Real_Vector (1 .. N);
582
      Diag_Adj  : Real_Vector (1 .. N);
583
 
584
      --  The vector Diag_Adj indicates the amount of change in each value,
585
      --  while Diag tracks the value itself and Values holds the values as
586
      --  they were at the beginning. As the changes typically will be small
587
      --  compared to the absolute value of Diag, at the end of each iteration
588
      --  Diag is computed as Diag + Diag_Adj thus avoiding accumulating
589
      --  rounding errors. This technique is due to Rutishauser.
590
 
591
   begin
592
      if Compute_Vectors
593
         and then (Vectors'Length (1) /= N or else Vectors'Length (2) /= N)
594
      then
595
         raise Constraint_Error with "incompatible matrix dimensions";
596
 
597
      elsif Values'Length /= N then
598
         raise Constraint_Error with "incompatible vector length";
599
 
600
      elsif not Is_Symmetric (M) then
601
         raise Constraint_Error with "matrix not symmetric";
602
      end if;
603
 
604
      --  Note: Only the locally declared matrix M and vectors (Diag, Diag_Adj)
605
      --        have lower bound equal to 1. The Vectors matrix may have
606
      --        different bounds, so take care indexing elements. Assignment
607
      --        as a whole is fine as sliding is automatic in that case.
608
 
609
      Vectors := (if not Compute_Vectors then (1 .. 0 => (1 .. 0 => 0.0))
610
                  else Unit_Matrix (Vectors'Length (1), Vectors'Length (2)));
611
      Values := Diagonal (M);
612
 
613
      Sweep : for Iteration in 1 .. Max_Iterations loop
614
 
615
         --  The first three iterations, perform rotation for any non-zero
616
         --  element. After this, rotate only for those that are not much
617
         --  smaller than the average off-diagnal element. After the fifth
618
         --  iteration, additionally zero out off-diagonal elements that are
619
         --  very small compared to elements on the diagonal with the same
620
         --  column or row index.
621
 
622
         Sum := Sum_Strict_Upper (M);
623
 
624
         exit Sweep when Sum = 0.0;
625
 
626
         Threshold := (if Iteration < 4 then 0.2 * Sum / Real (N**2) else 0.0);
627
 
628
         --  Iterate over all off-diagonal elements, rotating any that have
629
         --  an absolute value that exceeds the threshold.
630
 
631
         Diag := Values;
632
         Diag_Adj := (others => 0.0); -- Accumulates adjustments to Diag
633
 
634
         for Row in 1 .. N - 1 loop
635
            for Col in Row + 1 .. N loop
636
 
637
               --  If, before the rotation M (Row, Col) is tiny compared to
638
               --  Diag (Row) and Diag (Col), rotation is skipped. This is
639
               --  meaningful, as it produces no larger error than would be
640
               --  produced anyhow if the rotation had been performed.
641
               --  Suppress this optimization in the first four sweeps, so
642
               --  that this procedure can be used for computing eigenvectors
643
               --  of perturbed diagonal matrices.
644
 
645
               if Iteration > 4
646
                  and then Is_Tiny (M (Row, Col), Compared_To => Diag (Row))
647
                  and then Is_Tiny (M (Row, Col), Compared_To => Diag (Col))
648
               then
649
                  M (Row, Col) := 0.0;
650
 
651
               elsif abs M (Row, Col) > Threshold then
652
                  Perform_Rotation : declare
653
                     Tan : constant Real := Compute_Tan (M (Row, Col),
654
                                               Diag (Col) - Diag (Row));
655
                     Cos : constant Real := 1.0 / Sqrt (1.0 + Tan**2);
656
                     Sin : constant Real := Tan * Cos;
657
                     Tau : constant Real := Sin / (1.0 + Cos);
658
                     Adj : constant Real := Tan * M (Row, Col);
659
 
660
                  begin
661
                     Diag_Adj (Row) := Diag_Adj (Row) - Adj;
662
                     Diag_Adj (Col) := Diag_Adj (Col) + Adj;
663
                     Diag (Row) := Diag (Row) - Adj;
664
                     Diag (Col) := Diag (Col) + Adj;
665
 
666
                     M (Row, Col) := 0.0;
667
 
668
                     for J in 1 .. Row - 1 loop        --  1 <= J < Row
669
                        Rotate (M (J, Row), M (J, Col), Sin, Tau);
670
                     end loop;
671
 
672
                     for J in Row + 1 .. Col - 1 loop  --  Row < J < Col
673
                        Rotate (M (Row, J), M (J, Col), Sin, Tau);
674
                     end loop;
675
 
676
                     for J in Col + 1 .. N loop        --  Col < J <= N
677
                        Rotate (M (Row, J), M (Col, J), Sin, Tau);
678
                     end loop;
679
 
680
                     for J in Vectors'Range (1) loop
681
                        Rotate (Vectors (J, Row - 1 + Vectors'First (2)),
682
                                Vectors (J, Col - 1 + Vectors'First (2)),
683
                                Sin, Tau);
684
                     end loop;
685
                  end Perform_Rotation;
686
               end if;
687
            end loop;
688
         end loop;
689
 
690
         Values := Values + Diag_Adj;
691
      end loop Sweep;
692
 
693
      --  All normal matrices with valid values should converge perfectly.
694
 
695
      if Sum /= 0.0 then
696
         raise Constraint_Error with "eigensystem solution does not converge";
697
      end if;
698
   end Jacobi;
699
 
700
   -----------
701
   -- Solve --
702
   -----------
703
 
704
   function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector
705
      renames Instantiations.Solve;
706
 
707
   function Solve (A, X : Real_Matrix) return Real_Matrix
708
      renames Instantiations.Solve;
709
 
710
   ----------------------
711
   -- Sort_Eigensystem --
712
   ----------------------
713
 
714
   procedure Sort_Eigensystem
715
     (Values  : in out Real_Vector;
716
      Vectors : in out Real_Matrix)
717
   is
718
      procedure Swap (Left, Right : Integer);
719
      --  Swap Values (Left) with Values (Right), and also swap the
720
      --  corresponding eigenvectors. Note that lowerbounds may differ.
721
 
722
      function Less (Left, Right : Integer) return Boolean is
723
        (Values (Left) > Values (Right));
724
      --  Sort by decreasing eigenvalue, see RM G.3.1 (76).
725
 
726
      procedure Sort is new Generic_Anonymous_Array_Sort (Integer);
727
      --  Sorts eigenvalues and eigenvectors by decreasing value
728
 
729
      procedure Swap (Left, Right : Integer) is
730
      begin
731
         Swap (Values (Left), Values (Right));
732
         Swap_Column (Vectors, Left - Values'First + Vectors'First (2),
733
                               Right - Values'First + Vectors'First (2));
734
      end Swap;
735
 
736
   begin
737
      Sort (Values'First, Values'Last);
738
   end Sort_Eigensystem;
739
 
740
   ---------------
741
   -- Transpose --
742
   ---------------
743
 
744
   function Transpose (X : Real_Matrix) return Real_Matrix is
745
      R : Real_Matrix (X'Range (2), X'Range (1));
746
   begin
747
      Transpose (X, R);
748
      return R;
749
   end Transpose;
750
 
751
   -----------------
752
   -- Unit_Matrix --
753
   -----------------
754
 
755
   function Unit_Matrix
756
     (Order   : Positive;
757
      First_1 : Integer := 1;
758
      First_2 : Integer := 1) return Real_Matrix
759
     renames Instantiations.Unit_Matrix;
760
 
761
   -----------------
762
   -- Unit_Vector --
763
   -----------------
764
 
765
   function Unit_Vector
766
     (Index : Integer;
767
      Order : Positive;
768
      First : Integer := 1) return Real_Vector
769
     renames Instantiations.Unit_Vector;
770
 
771
end Ada.Numerics.Generic_Real_Arrays;

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.