|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT RUN-TIME COMPONENTS --
|
|
|
4 -- --
|
|
|
5 -- ADA.NUMERICS.GENERIC_REAL_ARRAYS --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
145
|
9 -- Copyright (C) 2006-2019, Free Software Foundation, Inc. --
|
|
111
|
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;
|
