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