?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fixing pr42337
[official-gcc.git] / gcc / ada / a-ngcoar.adb
blob9979d6baec6bba7413dc7cf73c255086dc4c77d7
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 2006-2009, 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 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
33 with System.Generic_Complex_BLAS;
34 with System.Generic_Complex_LAPACK;
35
36 package body Ada.Numerics.Generic_Complex_Arrays is
37
38 -- Operations involving inner products use BLAS library implementations.
39 -- This allows larger matrices and vectors to be computed efficiently,
40 -- taking into account memory hierarchy issues and vector instructions
41 -- that vary widely between machines.
42
43 -- Operations that are defined in terms of operations on the type Real,
44 -- such as addition, subtraction and scaling, are computed in the canonical
45 -- way looping over all elements.
46
47 -- Operations for solving linear systems and computing determinant,
48 -- eigenvalues, eigensystem and inverse, are implemented using the
49 -- LAPACK library.
50
51 type BLAS_Real_Vector is array (Integer range <>) of Real;
52
53 package BLAS is new System.Generic_Complex_BLAS
54 (Real => Real,
55 Complex_Types => Complex_Types,
56 Complex_Vector => Complex_Vector,
57 Complex_Matrix => Complex_Matrix);
58
59 package LAPACK is new System.Generic_Complex_LAPACK
60 (Real => Real,
61 Real_Vector => BLAS_Real_Vector,
62 Complex_Types => Complex_Types,
63 Complex_Vector => Complex_Vector,
64 Complex_Matrix => Complex_Matrix);
65
66 subtype Real is Real_Arrays.Real;
67 -- Work around visibility bug ???
68
69 use BLAS, LAPACK;
70
71 -- Procedure versions of functions returning unconstrained values.
72 -- This allows for inlining the function wrapper.
73
74 procedure Eigenvalues
75 (A : Complex_Matrix;
76 Values : out Real_Vector);
77
78 procedure Inverse
79 (A : Complex_Matrix;
80 R : out Complex_Matrix);
81
82 procedure Solve
83 (A : Complex_Matrix;
84 X : Complex_Vector;
85 B : out Complex_Vector);
86
87 procedure Solve
88 (A : Complex_Matrix;
89 X : Complex_Matrix;
90 B : out Complex_Matrix);
91
92 procedure Transpose is new System.Generic_Array_Operations.Transpose
93 (Scalar => Complex,
94 Matrix => Complex_Matrix);
95
96 -- Helper function that raises a Constraint_Error is the argument is
97 -- not a square matrix, and otherwise returns its length.
98
99 function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
100
101 -- Instantiating the following subprograms directly would lead to
102 -- name clashes, so use a local package.
103
104 package Instantiations is
105
106 ---------
107 -- "*" --
108 ---------
109
110 function "*" is new Vector_Scalar_Elementwise_Operation
111 (Left_Scalar => Complex,
112 Right_Scalar => Complex,
113 Result_Scalar => Complex,
114 Left_Vector => Complex_Vector,
115 Result_Vector => Complex_Vector,
116 Operation => "*");
117
118 function "*" is new Vector_Scalar_Elementwise_Operation
119 (Left_Scalar => Complex,
120 Right_Scalar => Real'Base,
121 Result_Scalar => Complex,
122 Left_Vector => Complex_Vector,
123 Result_Vector => Complex_Vector,
124 Operation => "*");
125
126 function "*" is new Scalar_Vector_Elementwise_Operation
127 (Left_Scalar => Complex,
128 Right_Scalar => Complex,
129 Result_Scalar => Complex,
130 Right_Vector => Complex_Vector,
131 Result_Vector => Complex_Vector,
132 Operation => "*");
133
134 function "*" is new Scalar_Vector_Elementwise_Operation
135 (Left_Scalar => Real'Base,
136 Right_Scalar => Complex,
137 Result_Scalar => Complex,
138 Right_Vector => Complex_Vector,
139 Result_Vector => Complex_Vector,
140 Operation => "*");
141
142 function "*" is new Inner_Product
143 (Left_Scalar => Complex,
144 Right_Scalar => Real'Base,
145 Result_Scalar => Complex,
146 Left_Vector => Complex_Vector,
147 Right_Vector => Real_Vector,
148 Zero => (0.0, 0.0));
149
150 function "*" is new Inner_Product
151 (Left_Scalar => Real'Base,
152 Right_Scalar => Complex,
153 Result_Scalar => Complex,
154 Left_Vector => Real_Vector,
155 Right_Vector => Complex_Vector,
156 Zero => (0.0, 0.0));
157
158 function "*" is new Outer_Product
159 (Left_Scalar => Complex,
160 Right_Scalar => Complex,
161 Result_Scalar => Complex,
162 Left_Vector => Complex_Vector,
163 Right_Vector => Complex_Vector,
164 Matrix => Complex_Matrix);
165
166 function "*" is new Outer_Product
167 (Left_Scalar => Real'Base,
168 Right_Scalar => Complex,
169 Result_Scalar => Complex,
170 Left_Vector => Real_Vector,
171 Right_Vector => Complex_Vector,
172 Matrix => Complex_Matrix);
173
174 function "*" is new Outer_Product
175 (Left_Scalar => Complex,
176 Right_Scalar => Real'Base,
177 Result_Scalar => Complex,
178 Left_Vector => Complex_Vector,
179 Right_Vector => Real_Vector,
180 Matrix => Complex_Matrix);
181
182 function "*" is new Matrix_Scalar_Elementwise_Operation
183 (Left_Scalar => Complex,
184 Right_Scalar => Complex,
185 Result_Scalar => Complex,
186 Left_Matrix => Complex_Matrix,
187 Result_Matrix => Complex_Matrix,
188 Operation => "*");
189
190 function "*" is new Matrix_Scalar_Elementwise_Operation
191 (Left_Scalar => Complex,
192 Right_Scalar => Real'Base,
193 Result_Scalar => Complex,
194 Left_Matrix => Complex_Matrix,
195 Result_Matrix => Complex_Matrix,
196 Operation => "*");
197
198 function "*" is new Scalar_Matrix_Elementwise_Operation
199 (Left_Scalar => Complex,
200 Right_Scalar => Complex,
201 Result_Scalar => Complex,
202 Right_Matrix => Complex_Matrix,
203 Result_Matrix => Complex_Matrix,
204 Operation => "*");
205
206 function "*" is new Scalar_Matrix_Elementwise_Operation
207 (Left_Scalar => Real'Base,
208 Right_Scalar => Complex,
209 Result_Scalar => Complex,
210 Right_Matrix => Complex_Matrix,
211 Result_Matrix => Complex_Matrix,
212 Operation => "*");
213
214 function "*" is new Matrix_Vector_Product
215 (Left_Scalar => Real'Base,
216 Right_Scalar => Complex,
217 Result_Scalar => Complex,
218 Matrix => Real_Matrix,
219 Right_Vector => Complex_Vector,
220 Result_Vector => Complex_Vector,
221 Zero => (0.0, 0.0));
222
223 function "*" is new Matrix_Vector_Product
224 (Left_Scalar => Complex,
225 Right_Scalar => Real'Base,
226 Result_Scalar => Complex,
227 Matrix => Complex_Matrix,
228 Right_Vector => Real_Vector,
229 Result_Vector => Complex_Vector,
230 Zero => (0.0, 0.0));
231
232 function "*" is new Vector_Matrix_Product
233 (Left_Scalar => Real'Base,
234 Right_Scalar => Complex,
235 Result_Scalar => Complex,
236 Left_Vector => Real_Vector,
237 Matrix => Complex_Matrix,
238 Result_Vector => Complex_Vector,
239 Zero => (0.0, 0.0));
240
241 function "*" is new Vector_Matrix_Product
242 (Left_Scalar => Complex,
243 Right_Scalar => Real'Base,
244 Result_Scalar => Complex,
245 Left_Vector => Complex_Vector,
246 Matrix => Real_Matrix,
247 Result_Vector => Complex_Vector,
248 Zero => (0.0, 0.0));
249
250 function "*" is new Matrix_Matrix_Product
251 (Left_Scalar => Real'Base,
252 Right_Scalar => Complex,
253 Result_Scalar => Complex,
254 Left_Matrix => Real_Matrix,
255 Right_Matrix => Complex_Matrix,
256 Result_Matrix => Complex_Matrix,
257 Zero => (0.0, 0.0));
258
259 function "*" is new Matrix_Matrix_Product
260 (Left_Scalar => Complex,
261 Right_Scalar => Real'Base,
262 Result_Scalar => Complex,
263 Left_Matrix => Complex_Matrix,
264 Right_Matrix => Real_Matrix,
265 Result_Matrix => Complex_Matrix,
266 Zero => (0.0, 0.0));
267
268 ---------
269 -- "+" --
270 ---------
271
272 function "+" is new Vector_Elementwise_Operation
273 (X_Scalar => Complex,
274 Result_Scalar => Complex,
275 X_Vector => Complex_Vector,
276 Result_Vector => Complex_Vector,
277 Operation => "+");
278
279 function "+" is new Vector_Vector_Elementwise_Operation
280 (Left_Scalar => Complex,
281 Right_Scalar => Complex,
282 Result_Scalar => Complex,
283 Left_Vector => Complex_Vector,
284 Right_Vector => Complex_Vector,
285 Result_Vector => Complex_Vector,
286 Operation => "+");
287
288 function "+" is new Vector_Vector_Elementwise_Operation
289 (Left_Scalar => Real'Base,
290 Right_Scalar => Complex,
291 Result_Scalar => Complex,
292 Left_Vector => Real_Vector,
293 Right_Vector => Complex_Vector,
294 Result_Vector => Complex_Vector,
295 Operation => "+");
296
297 function "+" is new Vector_Vector_Elementwise_Operation
298 (Left_Scalar => Complex,
299 Right_Scalar => Real'Base,
300 Result_Scalar => Complex,
301 Left_Vector => Complex_Vector,
302 Right_Vector => Real_Vector,
303 Result_Vector => Complex_Vector,
304 Operation => "+");
305
306 function "+" is new Matrix_Elementwise_Operation
307 (X_Scalar => Complex,
308 Result_Scalar => Complex,
309 X_Matrix => Complex_Matrix,
310 Result_Matrix => Complex_Matrix,
311 Operation => "+");
312
313 function "+" is new Matrix_Matrix_Elementwise_Operation
314 (Left_Scalar => Complex,
315 Right_Scalar => Complex,
316 Result_Scalar => Complex,
317 Left_Matrix => Complex_Matrix,
318 Right_Matrix => Complex_Matrix,
319 Result_Matrix => Complex_Matrix,
320 Operation => "+");
321
322 function "+" is new Matrix_Matrix_Elementwise_Operation
323 (Left_Scalar => Real'Base,
324 Right_Scalar => Complex,
325 Result_Scalar => Complex,
326 Left_Matrix => Real_Matrix,
327 Right_Matrix => Complex_Matrix,
328 Result_Matrix => Complex_Matrix,
329 Operation => "+");
330
331 function "+" is new Matrix_Matrix_Elementwise_Operation
332 (Left_Scalar => Complex,
333 Right_Scalar => Real'Base,
334 Result_Scalar => Complex,
335 Left_Matrix => Complex_Matrix,
336 Right_Matrix => Real_Matrix,
337 Result_Matrix => Complex_Matrix,
338 Operation => "+");
339
340 ---------
341 -- "-" --
342 ---------
343
344 function "-" is new Vector_Elementwise_Operation
345 (X_Scalar => Complex,
346 Result_Scalar => Complex,
347 X_Vector => Complex_Vector,
348 Result_Vector => Complex_Vector,
349 Operation => "-");
350
351 function "-" is new Vector_Vector_Elementwise_Operation
352 (Left_Scalar => Complex,
353 Right_Scalar => Complex,
354 Result_Scalar => Complex,
355 Left_Vector => Complex_Vector,
356 Right_Vector => Complex_Vector,
357 Result_Vector => Complex_Vector,
358 Operation => "-");
359
360 function "-" is new Vector_Vector_Elementwise_Operation
361 (Left_Scalar => Real'Base,
362 Right_Scalar => Complex,
363 Result_Scalar => Complex,
364 Left_Vector => Real_Vector,
365 Right_Vector => Complex_Vector,
366 Result_Vector => Complex_Vector,
367 Operation => "-");
368
369 function "-" is new Vector_Vector_Elementwise_Operation
370 (Left_Scalar => Complex,
371 Right_Scalar => Real'Base,
372 Result_Scalar => Complex,
373 Left_Vector => Complex_Vector,
374 Right_Vector => Real_Vector,
375 Result_Vector => Complex_Vector,
376 Operation => "-");
377
378 function "-" is new Matrix_Elementwise_Operation
379 (X_Scalar => Complex,
380 Result_Scalar => Complex,
381 X_Matrix => Complex_Matrix,
382 Result_Matrix => Complex_Matrix,
383 Operation => "-");
384
385 function "-" is new Matrix_Matrix_Elementwise_Operation
386 (Left_Scalar => Complex,
387 Right_Scalar => Complex,
388 Result_Scalar => Complex,
389 Left_Matrix => Complex_Matrix,
390 Right_Matrix => Complex_Matrix,
391 Result_Matrix => Complex_Matrix,
392 Operation => "-");
393
394 function "-" is new Matrix_Matrix_Elementwise_Operation
395 (Left_Scalar => Real'Base,
396 Right_Scalar => Complex,
397 Result_Scalar => Complex,
398 Left_Matrix => Real_Matrix,
399 Right_Matrix => Complex_Matrix,
400 Result_Matrix => Complex_Matrix,
401 Operation => "-");
402
403 function "-" is new Matrix_Matrix_Elementwise_Operation
404 (Left_Scalar => Complex,
405 Right_Scalar => Real'Base,
406 Result_Scalar => Complex,
407 Left_Matrix => Complex_Matrix,
408 Right_Matrix => Real_Matrix,
409 Result_Matrix => Complex_Matrix,
410 Operation => "-");
411
412 ---------
413 -- "/" --
414 ---------
415
416 function "/" is new Vector_Scalar_Elementwise_Operation
417 (Left_Scalar => Complex,
418 Right_Scalar => Complex,
419 Result_Scalar => Complex,
420 Left_Vector => Complex_Vector,
421 Result_Vector => Complex_Vector,
422 Operation => "/");
423
424 function "/" is new Vector_Scalar_Elementwise_Operation
425 (Left_Scalar => Complex,
426 Right_Scalar => Real'Base,
427 Result_Scalar => Complex,
428 Left_Vector => Complex_Vector,
429 Result_Vector => Complex_Vector,
430 Operation => "/");
431
432 function "/" is new Matrix_Scalar_Elementwise_Operation
433 (Left_Scalar => Complex,
434 Right_Scalar => Complex,
435 Result_Scalar => Complex,
436 Left_Matrix => Complex_Matrix,
437 Result_Matrix => Complex_Matrix,
438 Operation => "/");
439
440 function "/" is new Matrix_Scalar_Elementwise_Operation
441 (Left_Scalar => Complex,
442 Right_Scalar => Real'Base,
443 Result_Scalar => Complex,
444 Left_Matrix => Complex_Matrix,
445 Result_Matrix => Complex_Matrix,
446 Operation => "/");
447
448 --------------
449 -- Argument --
450 --------------
451
452 function Argument is new Vector_Elementwise_Operation
453 (X_Scalar => Complex,
454 Result_Scalar => Real'Base,
455 X_Vector => Complex_Vector,
456 Result_Vector => Real_Vector,
457 Operation => Argument);
458
459 function Argument is new Vector_Scalar_Elementwise_Operation
460 (Left_Scalar => Complex,
461 Right_Scalar => Real'Base,
462 Result_Scalar => Real'Base,
463 Left_Vector => Complex_Vector,
464 Result_Vector => Real_Vector,
465 Operation => Argument);
466
467 function Argument is new Matrix_Elementwise_Operation
468 (X_Scalar => Complex,
469 Result_Scalar => Real'Base,
470 X_Matrix => Complex_Matrix,
471 Result_Matrix => Real_Matrix,
472 Operation => Argument);
473
474 function Argument is new Matrix_Scalar_Elementwise_Operation
475 (Left_Scalar => Complex,
476 Right_Scalar => Real'Base,
477 Result_Scalar => Real'Base,
478 Left_Matrix => Complex_Matrix,
479 Result_Matrix => Real_Matrix,
480 Operation => Argument);
481
482 ----------------------------
483 -- Compose_From_Cartesian --
484 ----------------------------
485
486 function Compose_From_Cartesian is new Vector_Elementwise_Operation
487 (X_Scalar => Real'Base,
488 Result_Scalar => Complex,
489 X_Vector => Real_Vector,
490 Result_Vector => Complex_Vector,
491 Operation => Compose_From_Cartesian);
492
493 function Compose_From_Cartesian is
494 new Vector_Vector_Elementwise_Operation
495 (Left_Scalar => Real'Base,
496 Right_Scalar => Real'Base,
497 Result_Scalar => Complex,
498 Left_Vector => Real_Vector,
499 Right_Vector => Real_Vector,
500 Result_Vector => Complex_Vector,
501 Operation => Compose_From_Cartesian);
502
503 function Compose_From_Cartesian is new Matrix_Elementwise_Operation
504 (X_Scalar => Real'Base,
505 Result_Scalar => Complex,
506 X_Matrix => Real_Matrix,
507 Result_Matrix => Complex_Matrix,
508 Operation => Compose_From_Cartesian);
509
510 function Compose_From_Cartesian is
511 new Matrix_Matrix_Elementwise_Operation
512 (Left_Scalar => Real'Base,
513 Right_Scalar => Real'Base,
514 Result_Scalar => Complex,
515 Left_Matrix => Real_Matrix,
516 Right_Matrix => Real_Matrix,
517 Result_Matrix => Complex_Matrix,
518 Operation => Compose_From_Cartesian);
519
520 ------------------------
521 -- Compose_From_Polar --
522 ------------------------
523
524 function Compose_From_Polar is
525 new Vector_Vector_Elementwise_Operation
526 (Left_Scalar => Real'Base,
527 Right_Scalar => Real'Base,
528 Result_Scalar => Complex,
529 Left_Vector => Real_Vector,
530 Right_Vector => Real_Vector,
531 Result_Vector => Complex_Vector,
532 Operation => Compose_From_Polar);
533
534 function Compose_From_Polar is
535 new Vector_Vector_Scalar_Elementwise_Operation
536 (X_Scalar => Real'Base,
537 Y_Scalar => Real'Base,
538 Z_Scalar => Real'Base,
539 Result_Scalar => Complex,
540 X_Vector => Real_Vector,
541 Y_Vector => Real_Vector,
542 Result_Vector => Complex_Vector,
543 Operation => Compose_From_Polar);
544
545 function Compose_From_Polar is
546 new Matrix_Matrix_Elementwise_Operation
547 (Left_Scalar => Real'Base,
548 Right_Scalar => Real'Base,
549 Result_Scalar => Complex,
550 Left_Matrix => Real_Matrix,
551 Right_Matrix => Real_Matrix,
552 Result_Matrix => Complex_Matrix,
553 Operation => Compose_From_Polar);
554
555 function Compose_From_Polar is
556 new Matrix_Matrix_Scalar_Elementwise_Operation
557 (X_Scalar => Real'Base,
558 Y_Scalar => Real'Base,
559 Z_Scalar => Real'Base,
560 Result_Scalar => Complex,
561 X_Matrix => Real_Matrix,
562 Y_Matrix => Real_Matrix,
563 Result_Matrix => Complex_Matrix,
564 Operation => Compose_From_Polar);
565
566 ---------------
567 -- Conjugate --
568 ---------------
569
570 function Conjugate is new Vector_Elementwise_Operation
571 (X_Scalar => Complex,
572 Result_Scalar => Complex,
573 X_Vector => Complex_Vector,
574 Result_Vector => Complex_Vector,
575 Operation => Conjugate);
576
577 function Conjugate is new Matrix_Elementwise_Operation
578 (X_Scalar => Complex,
579 Result_Scalar => Complex,
580 X_Matrix => Complex_Matrix,
581 Result_Matrix => Complex_Matrix,
582 Operation => Conjugate);
583
584 --------
585 -- Im --
586 --------
587
588 function Im is new Vector_Elementwise_Operation
589 (X_Scalar => Complex,
590 Result_Scalar => Real'Base,
591 X_Vector => Complex_Vector,
592 Result_Vector => Real_Vector,
593 Operation => Im);
594
595 function Im is new Matrix_Elementwise_Operation
596 (X_Scalar => Complex,
597 Result_Scalar => Real'Base,
598 X_Matrix => Complex_Matrix,
599 Result_Matrix => Real_Matrix,
600 Operation => Im);
601
602 -------------
603 -- Modulus --
604 -------------
605
606 function Modulus is new Vector_Elementwise_Operation
607 (X_Scalar => Complex,
608 Result_Scalar => Real'Base,
609 X_Vector => Complex_Vector,
610 Result_Vector => Real_Vector,
611 Operation => Modulus);
612
613 function Modulus is new Matrix_Elementwise_Operation
614 (X_Scalar => Complex,
615 Result_Scalar => Real'Base,
616 X_Matrix => Complex_Matrix,
617 Result_Matrix => Real_Matrix,
618 Operation => Modulus);
619
620 --------
621 -- Re --
622 --------
623
624 function Re is new Vector_Elementwise_Operation
625 (X_Scalar => Complex,
626 Result_Scalar => Real'Base,
627 X_Vector => Complex_Vector,
628 Result_Vector => Real_Vector,
629 Operation => Re);
630
631 function Re is new Matrix_Elementwise_Operation
632 (X_Scalar => Complex,
633 Result_Scalar => Real'Base,
634 X_Matrix => Complex_Matrix,
635 Result_Matrix => Real_Matrix,
636 Operation => Re);
637
638 ------------
639 -- Set_Im --
640 ------------
641
642 procedure Set_Im is new Update_Vector_With_Vector
643 (X_Scalar => Complex,
644 Y_Scalar => Real'Base,
645 X_Vector => Complex_Vector,
646 Y_Vector => Real_Vector,
647 Update => Set_Im);
648
649 procedure Set_Im is new Update_Matrix_With_Matrix
650 (X_Scalar => Complex,
651 Y_Scalar => Real'Base,
652 X_Matrix => Complex_Matrix,
653 Y_Matrix => Real_Matrix,
654 Update => Set_Im);
655
656 ------------
657 -- Set_Re --
658 ------------
659
660 procedure Set_Re is new Update_Vector_With_Vector
661 (X_Scalar => Complex,
662 Y_Scalar => Real'Base,
663 X_Vector => Complex_Vector,
664 Y_Vector => Real_Vector,
665 Update => Set_Re);
666
667 procedure Set_Re is new Update_Matrix_With_Matrix
668 (X_Scalar => Complex,
669 Y_Scalar => Real'Base,
670 X_Matrix => Complex_Matrix,
671 Y_Matrix => Real_Matrix,
672 Update => Set_Re);
673
674 -----------------
675 -- Unit_Matrix --
676 -----------------
677
678 function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix
679 (Scalar => Complex,
680 Matrix => Complex_Matrix,
681 Zero => (0.0, 0.0),
682 One => (1.0, 0.0));
683
684 function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector
685 (Scalar => Complex,
686 Vector => Complex_Vector,
687 Zero => (0.0, 0.0),
688 One => (1.0, 0.0));
689
690 end Instantiations;
691
692 ---------
693 -- "*" --
694 ---------
695
696 function "*"
697 (Left : Complex_Vector;
698 Right : Complex_Vector) return Complex
699 is
700 begin
701 if Left'Length /= Right'Length then
702 raise Constraint_Error with
703 "vectors are of different length in inner product";
704 end if;
705
706 return dot (Left'Length, X => Left, Y => Right);
707 end "*";
708
709 function "*"
710 (Left : Real_Vector;
711 Right : Complex_Vector) return Complex
712 renames Instantiations."*";
713
714 function "*"
715 (Left : Complex_Vector;
716 Right : Real_Vector) return Complex
717 renames Instantiations."*";
718
719 function "*"
720 (Left : Complex;
721 Right : Complex_Vector) return Complex_Vector
722 renames Instantiations."*";
723
724 function "*"
725 (Left : Complex_Vector;
726 Right : Complex) return Complex_Vector
727 renames Instantiations."*";
728
729 function "*"
730 (Left : Real'Base;
731 Right : Complex_Vector) return Complex_Vector
732 renames Instantiations."*";
733
734 function "*"
735 (Left : Complex_Vector;
736 Right : Real'Base) return Complex_Vector
737 renames Instantiations."*";
738
739 function "*"
740 (Left : Complex_Matrix;
741 Right : Complex_Matrix)
742 return Complex_Matrix
743 is
744 R : Complex_Matrix (Left'Range (1), Right'Range (2));
745
746 begin
747 if Left'Length (2) /= Right'Length (1) then
748 raise Constraint_Error with
749 "incompatible dimensions in matrix-matrix multiplication";
750 end if;
751
752 gemm (Trans_A => No_Trans'Access,
753 Trans_B => No_Trans'Access,
754 M => Right'Length (2),
755 N => Left'Length (1),
756 K => Right'Length (1),
757 A => Right,
758 Ld_A => Right'Length (2),
759 B => Left,
760 Ld_B => Left'Length (2),
761 C => R,
762 Ld_C => R'Length (2));
763
764 return R;
765 end "*";
766
767 function "*"
768 (Left : Complex_Vector;
769 Right : Complex_Vector) return Complex_Matrix
770 renames Instantiations."*";
771
772 function "*"
773 (Left : Complex_Vector;
774 Right : Complex_Matrix) return Complex_Vector
775 is
776 R : Complex_Vector (Right'Range (2));
777
778 begin
779 if Left'Length /= Right'Length (1) then
780 raise Constraint_Error with
781 "incompatible dimensions in vector-matrix multiplication";
782 end if;
783
784 gemv (Trans => No_Trans'Access,
785 M => Right'Length (2),
786 N => Right'Length (1),
787 A => Right,
788 Ld_A => Right'Length (2),
789 X => Left,
790 Y => R);
791
792 return R;
793 end "*";
794
795 function "*"
796 (Left : Complex_Matrix;
797 Right : Complex_Vector) return Complex_Vector
798 is
799 R : Complex_Vector (Left'Range (1));
800
801 begin
802 if Left'Length (2) /= Right'Length then
803 raise Constraint_Error with
804 "incompatible dimensions in matrix-vector multiplication";
805 end if;
806
807 gemv (Trans => Trans'Access,
808 M => Left'Length (2),
809 N => Left'Length (1),
810 A => Left,
811 Ld_A => Left'Length (2),
812 X => Right,
813 Y => R);
814
815 return R;
816 end "*";
817
818 function "*"
819 (Left : Real_Matrix;
820 Right : Complex_Matrix) return Complex_Matrix
821 renames Instantiations."*";
822
823 function "*"
824 (Left : Complex_Matrix;
825 Right : Real_Matrix) return Complex_Matrix
826 renames Instantiations."*";
827
828 function "*"
829 (Left : Real_Vector;
830 Right : Complex_Vector) return Complex_Matrix
831 renames Instantiations."*";
832
833 function "*"
834 (Left : Complex_Vector;
835 Right : Real_Vector) return Complex_Matrix
836 renames Instantiations."*";
837
838 function "*"
839 (Left : Real_Vector;
840 Right : Complex_Matrix) return Complex_Vector
841 renames Instantiations."*";
842
843 function "*"
844 (Left : Complex_Vector;
845 Right : Real_Matrix) return Complex_Vector
846 renames Instantiations."*";
847
848 function "*"
849 (Left : Real_Matrix;
850 Right : Complex_Vector) return Complex_Vector
851 renames Instantiations."*";
852
853 function "*"
854 (Left : Complex_Matrix;
855 Right : Real_Vector) return Complex_Vector
856 renames Instantiations."*";
857
858 function "*"
859 (Left : Complex;
860 Right : Complex_Matrix) return Complex_Matrix
861 renames Instantiations."*";
862
863 function "*"
864 (Left : Complex_Matrix;
865 Right : Complex) return Complex_Matrix
866 renames Instantiations."*";
867
868 function "*"
869 (Left : Real'Base;
870 Right : Complex_Matrix) return Complex_Matrix
871 renames Instantiations."*";
872
873 function "*"
874 (Left : Complex_Matrix;
875 Right : Real'Base) return Complex_Matrix
876 renames Instantiations."*";
877
878 ---------
879 -- "+" --
880 ---------
881
882 function "+" (Right : Complex_Vector) return Complex_Vector
883 renames Instantiations."+";
884
885 function "+"
886 (Left : Complex_Vector;
887 Right : Complex_Vector) return Complex_Vector
888 renames Instantiations."+";
889
890 function "+"
891 (Left : Real_Vector;
892 Right : Complex_Vector) return Complex_Vector
893 renames Instantiations."+";
894
895 function "+"
896 (Left : Complex_Vector;
897 Right : Real_Vector) return Complex_Vector
898 renames Instantiations."+";
899
900 function "+" (Right : Complex_Matrix) return Complex_Matrix
901 renames Instantiations."+";
902
903 function "+"
904 (Left : Complex_Matrix;
905 Right : Complex_Matrix) return Complex_Matrix
906 renames Instantiations."+";
907
908 function "+"
909 (Left : Real_Matrix;
910 Right : Complex_Matrix) return Complex_Matrix
911 renames Instantiations."+";
912
913 function "+"
914 (Left : Complex_Matrix;
915 Right : Real_Matrix) return Complex_Matrix
916 renames Instantiations."+";
917
918 ---------
919 -- "-" --
920 ---------
921
922 function "-"
923 (Right : Complex_Vector) return Complex_Vector
924 renames Instantiations."-";
925
926 function "-"
927 (Left : Complex_Vector;
928 Right : Complex_Vector) return Complex_Vector
929 renames Instantiations."-";
930
931 function "-"
932 (Left : Real_Vector;
933 Right : Complex_Vector) return Complex_Vector
934 renames Instantiations."-";
935
936 function "-"
937 (Left : Complex_Vector;
938 Right : Real_Vector) return Complex_Vector
939 renames Instantiations."-";
940
941 function "-" (Right : Complex_Matrix) return Complex_Matrix
942 renames Instantiations."-";
943
944 function "-"
945 (Left : Complex_Matrix;
946 Right : Complex_Matrix) return Complex_Matrix
947 renames Instantiations."-";
948
949 function "-"
950 (Left : Real_Matrix;
951 Right : Complex_Matrix) return Complex_Matrix
952 renames Instantiations."-";
953
954 function "-"
955 (Left : Complex_Matrix;
956 Right : Real_Matrix) return Complex_Matrix
957 renames Instantiations."-";
958
959 ---------
960 -- "/" --
961 ---------
962
963 function "/"
964 (Left : Complex_Vector;
965 Right : Complex) return Complex_Vector
966 renames Instantiations."/";
967
968 function "/"
969 (Left : Complex_Vector;
970 Right : Real'Base) return Complex_Vector
971 renames Instantiations."/";
972
973 function "/"
974 (Left : Complex_Matrix;
975 Right : Complex) return Complex_Matrix
976 renames Instantiations."/";
977
978 function "/"
979 (Left : Complex_Matrix;
980 Right : Real'Base) return Complex_Matrix
981 renames Instantiations."/";
982
983 -----------
984 -- "abs" --
985 -----------
986
987 function "abs" (Right : Complex_Vector) return Complex is
988 begin
989 return (nrm2 (Right'Length, Right), 0.0);
990 end "abs";
991
992 --------------
993 -- Argument --
994 --------------
995
996 function Argument (X : Complex_Vector) return Real_Vector
997 renames Instantiations.Argument;
998
999 function Argument
1000 (X : Complex_Vector;
1001 Cycle : Real'Base) return Real_Vector
1002 renames Instantiations.Argument;
1003
1004 function Argument (X : Complex_Matrix) return Real_Matrix
1005 renames Instantiations.Argument;
1006
1007 function Argument
1008 (X : Complex_Matrix;
1009 Cycle : Real'Base) return Real_Matrix
1010 renames Instantiations.Argument;
1011
1012 ----------------------------
1013 -- Compose_From_Cartesian --
1014 ----------------------------
1015
1016 function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector
1017 renames Instantiations.Compose_From_Cartesian;
1018
1019 function Compose_From_Cartesian
1020 (Re : Real_Vector;
1021 Im : Real_Vector) return Complex_Vector
1022 renames Instantiations.Compose_From_Cartesian;
1023
1024 function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix
1025 renames Instantiations.Compose_From_Cartesian;
1026
1027 function Compose_From_Cartesian
1028 (Re : Real_Matrix;
1029 Im : Real_Matrix) return Complex_Matrix
1030 renames Instantiations.Compose_From_Cartesian;
1031
1032 ------------------------
1033 -- Compose_From_Polar --
1034 ------------------------
1035
1036 function Compose_From_Polar
1037 (Modulus : Real_Vector;
1038 Argument : Real_Vector) return Complex_Vector
1039 renames Instantiations.Compose_From_Polar;
1040
1041 function Compose_From_Polar
1042 (Modulus : Real_Vector;
1043 Argument : Real_Vector;
1044 Cycle : Real'Base) return Complex_Vector
1045 renames Instantiations.Compose_From_Polar;
1046
1047 function Compose_From_Polar
1048 (Modulus : Real_Matrix;
1049 Argument : Real_Matrix) return Complex_Matrix
1050 renames Instantiations.Compose_From_Polar;
1051
1052 function Compose_From_Polar
1053 (Modulus : Real_Matrix;
1054 Argument : Real_Matrix;
1055 Cycle : Real'Base) return Complex_Matrix
1056 renames Instantiations.Compose_From_Polar;
1057
1058 ---------------
1059 -- Conjugate --
1060 ---------------
1061
1062 function Conjugate (X : Complex_Vector) return Complex_Vector
1063 renames Instantiations.Conjugate;
1064
1065 function Conjugate (X : Complex_Matrix) return Complex_Matrix
1066 renames Instantiations.Conjugate;
1067
1068 -----------------
1069 -- Determinant --
1070 -----------------
1071
1072 function Determinant (A : Complex_Matrix) return Complex is
1073 N : constant Integer := Length (A);
1074 LU : Complex_Matrix (1 .. N, 1 .. N) := A;
1075 Piv : Integer_Vector (1 .. N);
1076 Info : aliased Integer := -1;
1077 Neg : Boolean;
1078 Det : Complex;
1079
1080 begin
1081 if N = 0 then
1082 return (0.0, 0.0);
1083 end if;
1084
1085 getrf (N, N, LU, N, Piv, Info'Access);
1086
1087 if Info /= 0 then
1088 raise Constraint_Error with "ill-conditioned matrix";
1089 end if;
1090
1091 Det := LU (1, 1);
1092 Neg := Piv (1) /= 1;
1093
1094 for J in 2 .. N loop
1095 Det := Det * LU (J, J);
1096 Neg := Neg xor (Piv (J) /= J);
1097 end loop;
1098
1099 if Neg then
1100 return -Det;
1101
1102 else
1103 return Det;
1104 end if;
1105 end Determinant;
1106
1107 -----------------
1108 -- Eigensystem --
1109 -----------------
1110
1111 procedure Eigensystem
1112 (A : Complex_Matrix;
1113 Values : out Real_Vector;
1114 Vectors : out Complex_Matrix)
1115 is
1116 Job_Z : aliased Character := 'V';
1117 Rng : aliased Character := 'A';
1118 Uplo : aliased Character := 'U';
1119
1120 N : constant Natural := Length (A);
1121 W : BLAS_Real_Vector (Values'Range);
1122 M : Integer;
1123 B : Complex_Matrix (1 .. N, 1 .. N);
1124 L_Work : Complex_Vector (1 .. 1);
1125 LR_Work : BLAS_Real_Vector (1 .. 1);
1126 LI_Work : Integer_Vector (1 .. 1);
1127 I_Supp_Z : Integer_Vector (1 .. 2 * N);
1128 Info : aliased Integer;
1129
1130 begin
1131 if Values'Length /= N then
1132 raise Constraint_Error with "wrong length for output vector";
1133 end if;
1134
1135 if Vectors'First (1) /= A'First (1)
1136 or else Vectors'Last (1) /= A'Last (1)
1137 or else Vectors'First (2) /= A'First (2)
1138 or else Vectors'Last (2) /= A'Last (2)
1139 then
1140 raise Constraint_Error with "wrong dimensions for output matrix";
1141 end if;
1142
1143 if N = 0 then
1144 return;
1145 end if;
1146
1147 -- Check for hermitian matrix ???
1148 -- Copy only required triangle ???
1149
1150 B := A;
1151
1152 -- Find size of work area
1153
1154 heevr
1155 (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1156 M => M,
1157 W => W,
1158 Z => Vectors,
1159 Ld_Z => N,
1160 I_Supp_Z => I_Supp_Z,
1161 Work => L_Work,
1162 L_Work => -1,
1163 R_Work => LR_Work,
1164 LR_Work => -1,
1165 I_Work => LI_Work,
1166 LI_Work => -1,
1167 Info => Info'Access);
1168
1169 if Info /= 0 then
1170 raise Constraint_Error;
1171 end if;
1172
1173 declare
1174 Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1175 R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
1176 I_Work : Integer_Vector (1 .. LI_Work (1));
1177
1178 begin
1179 heevr
1180 (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1181 M => M,
1182 W => W,
1183 Z => Vectors,
1184 Ld_Z => N,
1185 I_Supp_Z => I_Supp_Z,
1186 Work => Work,
1187 L_Work => Work'Length,
1188 R_Work => R_Work,
1189 LR_Work => LR_Work'Length,
1190 I_Work => I_Work,
1191 LI_Work => LI_Work'Length,
1192 Info => Info'Access);
1193
1194 if Info /= 0 then
1195 raise Constraint_Error with "inverting non-Hermitian matrix";
1196 end if;
1197
1198 for J in Values'Range loop
1199 Values (J) := W (J);
1200 end loop;
1201 end;
1202 end Eigensystem;
1203
1204 -----------------
1205 -- Eigenvalues --
1206 -----------------
1207
1208 procedure Eigenvalues
1209 (A : Complex_Matrix;
1210 Values : out Real_Vector)
1211 is
1212 Job_Z : aliased Character := 'N';
1213 Rng : aliased Character := 'A';
1214 Uplo : aliased Character := 'U';
1215 N : constant Natural := Length (A);
1216 B : Complex_Matrix (1 .. N, 1 .. N) := A;
1217 Z : Complex_Matrix (1 .. 1, 1 .. 1);
1218 W : BLAS_Real_Vector (Values'Range);
1219 L_Work : Complex_Vector (1 .. 1);
1220 LR_Work : BLAS_Real_Vector (1 .. 1);
1221 LI_Work : Integer_Vector (1 .. 1);
1222 I_Supp_Z : Integer_Vector (1 .. 2 * N);
1223 M : Integer;
1224 Info : aliased Integer;
1225
1226 begin
1227 if Values'Length /= N then
1228 raise Constraint_Error with "wrong length for output vector";
1229 end if;
1230
1231 if N = 0 then
1232 return;
1233 end if;
1234
1235 -- Check for hermitian matrix ???
1236
1237 -- Find size of work area
1238
1239 heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1240 M => M,
1241 W => W,
1242 Z => Z,
1243 Ld_Z => 1,
1244 I_Supp_Z => I_Supp_Z,
1245 Work => L_Work, L_Work => -1,
1246 R_Work => LR_Work, LR_Work => -1,
1247 I_Work => LI_Work, LI_Work => -1,
1248 Info => Info'Access);
1249
1250 if Info /= 0 then
1251 raise Constraint_Error;
1252 end if;
1253
1254 declare
1255 Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1256 R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
1257 I_Work : Integer_Vector (1 .. LI_Work (1));
1258 begin
1259 heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1260 M => M,
1261 W => W,
1262 Z => Z,
1263 Ld_Z => 1,
1264 I_Supp_Z => I_Supp_Z,
1265 Work => Work, L_Work => Work'Length,
1266 R_Work => R_Work, LR_Work => R_Work'Length,
1267 I_Work => I_Work, LI_Work => I_Work'Length,
1268 Info => Info'Access);
1269
1270 if Info /= 0 then
1271 raise Constraint_Error with "inverting singular matrix";
1272 end if;
1273
1274 for J in Values'Range loop
1275 Values (J) := W (J);
1276 end loop;
1277 end;
1278 end Eigenvalues;
1279
1280 function Eigenvalues (A : Complex_Matrix) return Real_Vector is
1281 R : Real_Vector (A'Range (1));
1282 begin
1283 Eigenvalues (A, R);
1284 return R;
1285 end Eigenvalues;
1286
1287 --------
1288 -- Im --
1289 --------
1290
1291 function Im (X : Complex_Vector) return Real_Vector
1292 renames Instantiations.Im;
1293
1294 function Im (X : Complex_Matrix) return Real_Matrix
1295 renames Instantiations.Im;
1296
1297 -------------
1298 -- Inverse --
1299 -------------
1300
1301 procedure Inverse (A : Complex_Matrix; R : out Complex_Matrix) is
1302 N : constant Integer := Length (A);
1303 Piv : Integer_Vector (1 .. N);
1304 L_Work : Complex_Vector (1 .. 1);
1305 Info : aliased Integer := -1;
1306
1307 begin
1308 -- All computations are done using column-major order, but this works
1309 -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
1310
1311 R := A;
1312
1313 -- Compute LU decomposition
1314
1315 getrf (M => N,
1316 N => N,
1317 A => R,
1318 Ld_A => N,
1319 I_Piv => Piv,
1320 Info => Info'Access);
1321
1322 if Info /= 0 then
1323 raise Constraint_Error with "inverting singular matrix";
1324 end if;
1325
1326 -- Determine size of work area
1327
1328 getri (N => N,
1329 A => R,
1330 Ld_A => N,
1331 I_Piv => Piv,
1332 Work => L_Work,
1333 L_Work => -1,
1334 Info => Info'Access);
1335
1336 if Info /= 0 then
1337 raise Constraint_Error;
1338 end if;
1339
1340 declare
1341 Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1342
1343 begin
1344 -- Compute inverse from LU decomposition
1345
1346 getri (N => N,
1347 A => R,
1348 Ld_A => N,
1349 I_Piv => Piv,
1350 Work => Work,
1351 L_Work => Work'Length,
1352 Info => Info'Access);
1353
1354 if Info /= 0 then
1355 raise Constraint_Error with "inverting singular matrix";
1356 end if;
1357
1358 -- ??? Should iterate with gerfs, based on implementation advice
1359 end;
1360 end Inverse;
1361
1362 function Inverse (A : Complex_Matrix) return Complex_Matrix is
1363 R : Complex_Matrix (A'Range (2), A'Range (1));
1364 begin
1365 Inverse (A, R);
1366 return R;
1367 end Inverse;
1368
1369 -------------
1370 -- Modulus --
1371 -------------
1372
1373 function Modulus (X : Complex_Vector) return Real_Vector
1374 renames Instantiations.Modulus;
1375
1376 function Modulus (X : Complex_Matrix) return Real_Matrix
1377 renames Instantiations.Modulus;
1378
1379 --------
1380 -- Re --
1381 --------
1382
1383 function Re (X : Complex_Vector) return Real_Vector
1384 renames Instantiations.Re;
1385
1386 function Re (X : Complex_Matrix) return Real_Matrix
1387 renames Instantiations.Re;
1388
1389 ------------
1390 -- Set_Im --
1391 ------------
1392
1393 procedure Set_Im
1394 (X : in out Complex_Matrix;
1395 Im : Real_Matrix)
1396 renames Instantiations.Set_Im;
1397
1398 procedure Set_Im
1399 (X : in out Complex_Vector;
1400 Im : Real_Vector)
1401 renames Instantiations.Set_Im;
1402
1403 ------------
1404 -- Set_Re --
1405 ------------
1406
1407 procedure Set_Re
1408 (X : in out Complex_Matrix;
1409 Re : Real_Matrix)
1410 renames Instantiations.Set_Re;
1411
1412 procedure Set_Re
1413 (X : in out Complex_Vector;
1414 Re : Real_Vector)
1415 renames Instantiations.Set_Re;
1416
1417 -----------
1418 -- Solve --
1419 -----------
1420
1421 procedure Solve
1422 (A : Complex_Matrix;
1423 X : Complex_Vector;
1424 B : out Complex_Vector)
1425 is
1426 begin
1427 if Length (A) /= X'Length then
1428 raise Constraint_Error with
1429 "incompatible matrix and vector dimensions";
1430 end if;
1431
1432 -- ??? Should solve directly, is faster and more accurate
1433
1434 B := Inverse (A) * X;
1435 end Solve;
1436
1437 procedure Solve
1438 (A : Complex_Matrix;
1439 X : Complex_Matrix;
1440 B : out Complex_Matrix)
1441 is
1442 begin
1443 if Length (A) /= X'Length (1) then
1444 raise Constraint_Error with "incompatible matrix dimensions";
1445 end if;
1446
1447 -- ??? Should solve directly, is faster and more accurate
1448
1449 B := Inverse (A) * X;
1450 end Solve;
1451
1452 function Solve
1453 (A : Complex_Matrix;
1454 X : Complex_Vector) return Complex_Vector
1455 is
1456 B : Complex_Vector (A'Range (2));
1457 begin
1458 Solve (A, X, B);
1459 return B;
1460 end Solve;
1461
1462 function Solve (A, X : Complex_Matrix) return Complex_Matrix is
1463 B : Complex_Matrix (A'Range (2), X'Range (2));
1464 begin
1465 Solve (A, X, B);
1466 return B;
1467 end Solve;
1468
1469 ---------------
1470 -- Transpose --
1471 ---------------
1472
1473 function Transpose
1474 (X : Complex_Matrix) return Complex_Matrix
1475 is
1476 R : Complex_Matrix (X'Range (2), X'Range (1));
1477 begin
1478 Transpose (X, R);
1479 return R;
1480 end Transpose;
1481
1482 -----------------
1483 -- Unit_Matrix --
1484 -----------------
1485
1486 function Unit_Matrix
1487 (Order : Positive;
1488 First_1 : Integer := 1;
1489 First_2 : Integer := 1) return Complex_Matrix
1490 renames Instantiations.Unit_Matrix;
1491
1492 -----------------
1493 -- Unit_Vector --
1494 -----------------
1495
1496 function Unit_Vector
1497 (Index : Integer;
1498 Order : Positive;
1499 First : Integer := 1) return Complex_Vector
1500 renames Instantiations.Unit_Vector;
1501
1502 end Ada.Numerics.Generic_Complex_Arrays;
Mirror of the offical gcc git repository hosted at gcc.gnu.org