1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS --
9 -- Copyright (C) 2006-2009, Free Software Foundation, Inc. --
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. --
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. --
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/>. --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 ------------------------------------------------------------------------------
32 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
33 with System.Generic_Complex_BLAS;
34 with System.Generic_Complex_LAPACK;
36 package body Ada.Numerics.Generic_Complex_Arrays is
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.
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.
47 -- Operations for solving linear systems and computing determinant,
48 -- eigenvalues, eigensystem and inverse, are implemented using the
51 type BLAS_Real_Vector is array (Integer range <>) of Real;
53 package BLAS is new System.Generic_Complex_BLAS
55 Complex_Types => Complex_Types,
56 Complex_Vector => Complex_Vector,
57 Complex_Matrix => Complex_Matrix);
59 package LAPACK is new System.Generic_Complex_LAPACK
61 Real_Vector => BLAS_Real_Vector,
62 Complex_Types => Complex_Types,
63 Complex_Vector => Complex_Vector,
64 Complex_Matrix => Complex_Matrix);
66 subtype Real is Real_Arrays.Real;
67 -- Work around visibility bug ???
71 -- Procedure versions of functions returning unconstrained values.
72 -- This allows for inlining the function wrapper.
76 Values : out Real_Vector);
80 R : out Complex_Matrix);
85 B : out Complex_Vector);
90 B : out Complex_Matrix);
92 procedure Transpose is new System.Generic_Array_Operations.Transpose
94 Matrix => Complex_Matrix);
96 -- Helper function that raises a Constraint_Error is the argument is
97 -- not a square matrix, and otherwise returns its length.
99 function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
101 -- Instantiating the following subprograms directly would lead to
102 -- name clashes, so use a local package.
104 package Instantiations is
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,
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,
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,
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,
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,
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,
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);
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);
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);
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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);
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);
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);
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);
482 ----------------------------
483 -- Compose_From_Cartesian --
484 ----------------------------
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);
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);
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);
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);
520 ------------------------
521 -- Compose_From_Polar --
522 ------------------------
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);
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);
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);
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);
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);
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);
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,
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,
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);
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);
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,
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,
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,
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,
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,
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,
678 function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix
680 Matrix => Complex_Matrix,
684 function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector
686 Vector => Complex_Vector,
697 (Left : Complex_Vector;
698 Right : Complex_Vector) return Complex
701 if Left'Length /= Right'Length then
702 raise Constraint_Error with
703 "vectors are of different length in inner product";
706 return dot (Left'Length, X => Left, Y => Right);
711 Right : Complex_Vector) return Complex
712 renames Instantiations."*";
715 (Left : Complex_Vector;
716 Right : Real_Vector) return Complex
717 renames Instantiations."*";
721 Right : Complex_Vector) return Complex_Vector
722 renames Instantiations."*";
725 (Left : Complex_Vector;
726 Right : Complex) return Complex_Vector
727 renames Instantiations."*";
731 Right : Complex_Vector) return Complex_Vector
732 renames Instantiations."*";
735 (Left : Complex_Vector;
736 Right : Real'Base) return Complex_Vector
737 renames Instantiations."*";
740 (Left : Complex_Matrix;
741 Right : Complex_Matrix)
742 return Complex_Matrix
744 R : Complex_Matrix (Left'Range (1), Right'Range (2));
747 if Left'Length (2) /= Right'Length (1) then
748 raise Constraint_Error with
749 "incompatible dimensions in matrix-matrix multiplication";
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),
758 Ld_A => Right'Length (2),
760 Ld_B => Left'Length (2),
762 Ld_C => R'Length (2));
768 (Left : Complex_Vector;
769 Right : Complex_Vector) return Complex_Matrix
770 renames Instantiations."*";
773 (Left : Complex_Vector;
774 Right : Complex_Matrix) return Complex_Vector
776 R : Complex_Vector (Right'Range (2));
779 if Left'Length /= Right'Length (1) then
780 raise Constraint_Error with
781 "incompatible dimensions in vector-matrix multiplication";
784 gemv (Trans => No_Trans'Access,
785 M => Right'Length (2),
786 N => Right'Length (1),
788 Ld_A => Right'Length (2),
796 (Left : Complex_Matrix;
797 Right : Complex_Vector) return Complex_Vector
799 R : Complex_Vector (Left'Range (1));
802 if Left'Length (2) /= Right'Length then
803 raise Constraint_Error with
804 "incompatible dimensions in matrix-vector multiplication";
807 gemv (Trans => Trans'Access,
808 M => Left'Length (2),
809 N => Left'Length (1),
811 Ld_A => Left'Length (2),
820 Right : Complex_Matrix) return Complex_Matrix
821 renames Instantiations."*";
824 (Left : Complex_Matrix;
825 Right : Real_Matrix) return Complex_Matrix
826 renames Instantiations."*";
830 Right : Complex_Vector) return Complex_Matrix
831 renames Instantiations."*";
834 (Left : Complex_Vector;
835 Right : Real_Vector) return Complex_Matrix
836 renames Instantiations."*";
840 Right : Complex_Matrix) return Complex_Vector
841 renames Instantiations."*";
844 (Left : Complex_Vector;
845 Right : Real_Matrix) return Complex_Vector
846 renames Instantiations."*";
850 Right : Complex_Vector) return Complex_Vector
851 renames Instantiations."*";
854 (Left : Complex_Matrix;
855 Right : Real_Vector) return Complex_Vector
856 renames Instantiations."*";
860 Right : Complex_Matrix) return Complex_Matrix
861 renames Instantiations."*";
864 (Left : Complex_Matrix;
865 Right : Complex) return Complex_Matrix
866 renames Instantiations."*";
870 Right : Complex_Matrix) return Complex_Matrix
871 renames Instantiations."*";
874 (Left : Complex_Matrix;
875 Right : Real'Base) return Complex_Matrix
876 renames Instantiations."*";
882 function "+" (Right : Complex_Vector) return Complex_Vector
883 renames Instantiations."+";
886 (Left : Complex_Vector;
887 Right : Complex_Vector) return Complex_Vector
888 renames Instantiations."+";
892 Right : Complex_Vector) return Complex_Vector
893 renames Instantiations."+";
896 (Left : Complex_Vector;
897 Right : Real_Vector) return Complex_Vector
898 renames Instantiations."+";
900 function "+" (Right : Complex_Matrix) return Complex_Matrix
901 renames Instantiations."+";
904 (Left : Complex_Matrix;
905 Right : Complex_Matrix) return Complex_Matrix
906 renames Instantiations."+";
910 Right : Complex_Matrix) return Complex_Matrix
911 renames Instantiations."+";
914 (Left : Complex_Matrix;
915 Right : Real_Matrix) return Complex_Matrix
916 renames Instantiations."+";
923 (Right : Complex_Vector) return Complex_Vector
924 renames Instantiations."-";
927 (Left : Complex_Vector;
928 Right : Complex_Vector) return Complex_Vector
929 renames Instantiations."-";
933 Right : Complex_Vector) return Complex_Vector
934 renames Instantiations."-";
937 (Left : Complex_Vector;
938 Right : Real_Vector) return Complex_Vector
939 renames Instantiations."-";
941 function "-" (Right : Complex_Matrix) return Complex_Matrix
942 renames Instantiations."-";
945 (Left : Complex_Matrix;
946 Right : Complex_Matrix) return Complex_Matrix
947 renames Instantiations."-";
951 Right : Complex_Matrix) return Complex_Matrix
952 renames Instantiations."-";
955 (Left : Complex_Matrix;
956 Right : Real_Matrix) return Complex_Matrix
957 renames Instantiations."-";
964 (Left : Complex_Vector;
965 Right : Complex) return Complex_Vector
966 renames Instantiations."/";
969 (Left : Complex_Vector;
970 Right : Real'Base) return Complex_Vector
971 renames Instantiations."/";
974 (Left : Complex_Matrix;
975 Right : Complex) return Complex_Matrix
976 renames Instantiations."/";
979 (Left : Complex_Matrix;
980 Right : Real'Base) return Complex_Matrix
981 renames Instantiations."/";
987 function "abs" (Right : Complex_Vector) return Complex is
989 return (nrm2 (Right'Length, Right), 0.0);
996 function Argument (X : Complex_Vector) return Real_Vector
997 renames Instantiations.Argument;
1000 (X : Complex_Vector;
1001 Cycle : Real'Base) return Real_Vector
1002 renames Instantiations.Argument;
1004 function Argument (X : Complex_Matrix) return Real_Matrix
1005 renames Instantiations.Argument;
1008 (X : Complex_Matrix;
1009 Cycle : Real'Base) return Real_Matrix
1010 renames Instantiations.Argument;
1012 ----------------------------
1013 -- Compose_From_Cartesian --
1014 ----------------------------
1016 function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector
1017 renames Instantiations.Compose_From_Cartesian;
1019 function Compose_From_Cartesian
1021 Im : Real_Vector) return Complex_Vector
1022 renames Instantiations.Compose_From_Cartesian;
1024 function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix
1025 renames Instantiations.Compose_From_Cartesian;
1027 function Compose_From_Cartesian
1029 Im : Real_Matrix) return Complex_Matrix
1030 renames Instantiations.Compose_From_Cartesian;
1032 ------------------------
1033 -- Compose_From_Polar --
1034 ------------------------
1036 function Compose_From_Polar
1037 (Modulus : Real_Vector;
1038 Argument : Real_Vector) return Complex_Vector
1039 renames Instantiations.Compose_From_Polar;
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;
1047 function Compose_From_Polar
1048 (Modulus : Real_Matrix;
1049 Argument : Real_Matrix) return Complex_Matrix
1050 renames Instantiations.Compose_From_Polar;
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;
1062 function Conjugate (X : Complex_Vector) return Complex_Vector
1063 renames Instantiations.Conjugate;
1065 function Conjugate (X : Complex_Matrix) return Complex_Matrix
1066 renames Instantiations.Conjugate;
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;
1085 getrf (N, N, LU, N, Piv, Info'Access);
1088 raise Constraint_Error with "ill-conditioned matrix";
1092 Neg := Piv (1) /= 1;
1094 for J in 2 .. N loop
1095 Det := Det * LU (J, J);
1096 Neg := Neg xor (Piv (J) /= J);
1111 procedure Eigensystem
1112 (A : Complex_Matrix;
1113 Values : out Real_Vector;
1114 Vectors : out Complex_Matrix)
1116 Job_Z : aliased Character := 'V';
1117 Rng : aliased Character := 'A';
1118 Uplo : aliased Character := 'U';
1120 N : constant Natural := Length (A);
1121 W : BLAS_Real_Vector (Values'Range);
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;
1131 if Values'Length /= N then
1132 raise Constraint_Error with "wrong length for output vector";
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)
1140 raise Constraint_Error with "wrong dimensions for output matrix";
1147 -- Check for hermitian matrix ???
1148 -- Copy only required triangle ???
1152 -- Find size of work area
1155 (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1160 I_Supp_Z => I_Supp_Z,
1167 Info => Info'Access);
1170 raise Constraint_Error;
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));
1180 (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1185 I_Supp_Z => I_Supp_Z,
1187 L_Work => Work'Length,
1189 LR_Work => LR_Work'Length,
1191 LI_Work => LI_Work'Length,
1192 Info => Info'Access);
1195 raise Constraint_Error with "inverting non-Hermitian matrix";
1198 for J in Values'Range loop
1199 Values (J) := W (J);
1208 procedure Eigenvalues
1209 (A : Complex_Matrix;
1210 Values : out Real_Vector)
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);
1224 Info : aliased Integer;
1227 if Values'Length /= N then
1228 raise Constraint_Error with "wrong length for output vector";
1235 -- Check for hermitian matrix ???
1237 -- Find size of work area
1239 heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
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);
1251 raise Constraint_Error;
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));
1259 heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
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);
1271 raise Constraint_Error with "inverting singular matrix";
1274 for J in Values'Range loop
1275 Values (J) := W (J);
1280 function Eigenvalues (A : Complex_Matrix) return Real_Vector is
1281 R : Real_Vector (A'Range (1));
1291 function Im (X : Complex_Vector) return Real_Vector
1292 renames Instantiations.Im;
1294 function Im (X : Complex_Matrix) return Real_Matrix
1295 renames Instantiations.Im;
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;
1308 -- All computations are done using column-major order, but this works
1309 -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
1313 -- Compute LU decomposition
1320 Info => Info'Access);
1323 raise Constraint_Error with "inverting singular matrix";
1326 -- Determine size of work area
1334 Info => Info'Access);
1337 raise Constraint_Error;
1341 Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1344 -- Compute inverse from LU decomposition
1351 L_Work => Work'Length,
1352 Info => Info'Access);
1355 raise Constraint_Error with "inverting singular matrix";
1358 -- ??? Should iterate with gerfs, based on implementation advice
1362 function Inverse (A : Complex_Matrix) return Complex_Matrix is
1363 R : Complex_Matrix (A'Range (2), A'Range (1));
1373 function Modulus (X : Complex_Vector) return Real_Vector
1374 renames Instantiations.Modulus;
1376 function Modulus (X : Complex_Matrix) return Real_Matrix
1377 renames Instantiations.Modulus;
1383 function Re (X : Complex_Vector) return Real_Vector
1384 renames Instantiations.Re;
1386 function Re (X : Complex_Matrix) return Real_Matrix
1387 renames Instantiations.Re;
1394 (X : in out Complex_Matrix;
1396 renames Instantiations.Set_Im;
1399 (X : in out Complex_Vector;
1401 renames Instantiations.Set_Im;
1408 (X : in out Complex_Matrix;
1410 renames Instantiations.Set_Re;
1413 (X : in out Complex_Vector;
1415 renames Instantiations.Set_Re;
1422 (A : Complex_Matrix;
1424 B : out Complex_Vector)
1427 if Length (A) /= X'Length then
1428 raise Constraint_Error with
1429 "incompatible matrix and vector dimensions";
1432 -- ??? Should solve directly, is faster and more accurate
1434 B := Inverse (A) * X;
1438 (A : Complex_Matrix;
1440 B : out Complex_Matrix)
1443 if Length (A) /= X'Length (1) then
1444 raise Constraint_Error with "incompatible matrix dimensions";
1447 -- ??? Should solve directly, is faster and more accurate
1449 B := Inverse (A) * X;
1453 (A : Complex_Matrix;
1454 X : Complex_Vector) return Complex_Vector
1456 B : Complex_Vector (A'Range (2));
1462 function Solve (A, X : Complex_Matrix) return Complex_Matrix is
1463 B : Complex_Matrix (A'Range (2), X'Range (2));
1474 (X : Complex_Matrix) return Complex_Matrix
1476 R : Complex_Matrix (X'Range (2), X'Range (1));
1486 function Unit_Matrix
1488 First_1 : Integer := 1;
1489 First_2 : Integer := 1) return Complex_Matrix
1490 renames Instantiations.Unit_Matrix;
1496 function Unit_Vector
1499 First : Integer := 1) return Complex_Vector
1500 renames Instantiations.Unit_Vector;
1502 end Ada.Numerics.Generic_Complex_Arrays;