1 ------------------------------------------------------------------------------
3 -- GNAT COMPILER COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS --
9 -- Copyright (C) 2006-2007, 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 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. --
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. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
35 with System.Generic_Complex_BLAS;
36 with System.Generic_Complex_LAPACK;
38 package body Ada.Numerics.Generic_Complex_Arrays is
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.
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.
49 -- Operations for solving linear systems and computing determinant,
50 -- eigenvalues, eigensystem and inverse, are implemented using the
53 type BLAS_Real_Vector is array (Integer range <>) of Real;
55 package BLAS is new System.Generic_Complex_BLAS
57 Complex_Types => Complex_Types,
58 Complex_Vector => Complex_Vector,
59 Complex_Matrix => Complex_Matrix);
61 package LAPACK is new System.Generic_Complex_LAPACK
63 Real_Vector => BLAS_Real_Vector,
64 Complex_Types => Complex_Types,
65 Complex_Vector => Complex_Vector,
66 Complex_Matrix => Complex_Matrix);
68 subtype Real is Real_Arrays.Real;
69 -- Work around visibility bug ???
73 -- Procedure versions of functions returning unconstrained values.
74 -- This allows for inlining the function wrapper.
78 Values : out Real_Vector);
82 R : out Complex_Matrix);
87 B : out Complex_Vector);
92 B : out Complex_Matrix);
94 procedure Transpose is new System.Generic_Array_Operations.Transpose
96 Matrix => Complex_Matrix);
98 -- Helper function that raises a Constraint_Error is the argument is
99 -- not a square matrix, and otherwise returns its length.
101 function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
103 -- Instantiating the following subprograms directly would lead to
104 -- name clashes, so use a local package.
106 package Instantiations is
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,
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,
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,
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,
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,
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,
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);
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);
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);
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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);
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);
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);
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);
484 ----------------------------
485 -- Compose_From_Cartesian --
486 ----------------------------
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);
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);
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);
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);
522 ------------------------
523 -- Compose_From_Polar --
524 ------------------------
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);
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);
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);
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);
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);
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);
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,
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,
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);
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);
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,
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,
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,
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,
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,
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,
680 function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix
682 Matrix => Complex_Matrix,
686 function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector
688 Vector => Complex_Vector,
699 (Left : Complex_Vector;
700 Right : Complex_Vector) return Complex
703 if Left'Length /= Right'Length then
704 raise Constraint_Error with
705 "vectors are of different length in inner product";
708 return dot (Left'Length, X => Left, Y => Right);
713 Right : Complex_Vector) return Complex
714 renames Instantiations."*";
717 (Left : Complex_Vector;
718 Right : Real_Vector) return Complex
719 renames Instantiations."*";
723 Right : Complex_Vector) return Complex_Vector
724 renames Instantiations."*";
727 (Left : Complex_Vector;
728 Right : Complex) return Complex_Vector
729 renames Instantiations."*";
733 Right : Complex_Vector) return Complex_Vector
734 renames Instantiations."*";
737 (Left : Complex_Vector;
738 Right : Real'Base) return Complex_Vector
739 renames Instantiations."*";
742 (Left : Complex_Matrix;
743 Right : Complex_Matrix)
744 return Complex_Matrix
746 R : Complex_Matrix (Left'Range (1), Right'Range (2));
749 if Left'Length (2) /= Right'Length (1) then
750 raise Constraint_Error with
751 "incompatible dimensions in matrix-matrix multiplication";
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),
760 Ld_A => Right'Length (2),
762 Ld_B => Left'Length (2),
764 Ld_C => R'Length (2));
770 (Left : Complex_Vector;
771 Right : Complex_Vector) return Complex_Matrix
772 renames Instantiations."*";
775 (Left : Complex_Vector;
776 Right : Complex_Matrix) return Complex_Vector
778 R : Complex_Vector (Right'Range (2));
781 if Left'Length /= Right'Length (1) then
782 raise Constraint_Error with
783 "incompatible dimensions in vector-matrix multiplication";
786 gemv (Trans => No_Trans'Access,
787 M => Right'Length (2),
788 N => Right'Length (1),
790 Ld_A => Right'Length (2),
798 (Left : Complex_Matrix;
799 Right : Complex_Vector) return Complex_Vector
801 R : Complex_Vector (Left'Range (1));
804 if Left'Length (2) /= Right'Length then
805 raise Constraint_Error with
806 "incompatible dimensions in matrix-vector multiplication";
809 gemv (Trans => Trans'Access,
810 M => Left'Length (2),
811 N => Left'Length (1),
813 Ld_A => Left'Length (2),
822 Right : Complex_Matrix) return Complex_Matrix
823 renames Instantiations."*";
826 (Left : Complex_Matrix;
827 Right : Real_Matrix) return Complex_Matrix
828 renames Instantiations."*";
832 Right : Complex_Vector) return Complex_Matrix
833 renames Instantiations."*";
836 (Left : Complex_Vector;
837 Right : Real_Vector) return Complex_Matrix
838 renames Instantiations."*";
842 Right : Complex_Matrix) return Complex_Vector
843 renames Instantiations."*";
846 (Left : Complex_Vector;
847 Right : Real_Matrix) return Complex_Vector
848 renames Instantiations."*";
852 Right : Complex_Vector) return Complex_Vector
853 renames Instantiations."*";
856 (Left : Complex_Matrix;
857 Right : Real_Vector) return Complex_Vector
858 renames Instantiations."*";
862 Right : Complex_Matrix) return Complex_Matrix
863 renames Instantiations."*";
866 (Left : Complex_Matrix;
867 Right : Complex) return Complex_Matrix
868 renames Instantiations."*";
872 Right : Complex_Matrix) return Complex_Matrix
873 renames Instantiations."*";
876 (Left : Complex_Matrix;
877 Right : Real'Base) return Complex_Matrix
878 renames Instantiations."*";
884 function "+" (Right : Complex_Vector) return Complex_Vector
885 renames Instantiations."+";
888 (Left : Complex_Vector;
889 Right : Complex_Vector) return Complex_Vector
890 renames Instantiations."+";
894 Right : Complex_Vector) return Complex_Vector
895 renames Instantiations."+";
898 (Left : Complex_Vector;
899 Right : Real_Vector) return Complex_Vector
900 renames Instantiations."+";
902 function "+" (Right : Complex_Matrix) return Complex_Matrix
903 renames Instantiations."+";
906 (Left : Complex_Matrix;
907 Right : Complex_Matrix) return Complex_Matrix
908 renames Instantiations."+";
912 Right : Complex_Matrix) return Complex_Matrix
913 renames Instantiations."+";
916 (Left : Complex_Matrix;
917 Right : Real_Matrix) return Complex_Matrix
918 renames Instantiations."+";
925 (Right : Complex_Vector) return Complex_Vector
926 renames Instantiations."-";
929 (Left : Complex_Vector;
930 Right : Complex_Vector) return Complex_Vector
931 renames Instantiations."-";
935 Right : Complex_Vector) return Complex_Vector
936 renames Instantiations."-";
939 (Left : Complex_Vector;
940 Right : Real_Vector) return Complex_Vector
941 renames Instantiations."-";
943 function "-" (Right : Complex_Matrix) return Complex_Matrix
944 renames Instantiations."-";
947 (Left : Complex_Matrix;
948 Right : Complex_Matrix) return Complex_Matrix
949 renames Instantiations."-";
953 Right : Complex_Matrix) return Complex_Matrix
954 renames Instantiations."-";
957 (Left : Complex_Matrix;
958 Right : Real_Matrix) return Complex_Matrix
959 renames Instantiations."-";
966 (Left : Complex_Vector;
967 Right : Complex) return Complex_Vector
968 renames Instantiations."/";
971 (Left : Complex_Vector;
972 Right : Real'Base) return Complex_Vector
973 renames Instantiations."/";
976 (Left : Complex_Matrix;
977 Right : Complex) return Complex_Matrix
978 renames Instantiations."/";
981 (Left : Complex_Matrix;
982 Right : Real'Base) return Complex_Matrix
983 renames Instantiations."/";
989 function "abs" (Right : Complex_Vector) return Complex is
991 return (nrm2 (Right'Length, Right), 0.0);
998 function Argument (X : Complex_Vector) return Real_Vector
999 renames Instantiations.Argument;
1002 (X : Complex_Vector;
1003 Cycle : Real'Base) return Real_Vector
1004 renames Instantiations.Argument;
1006 function Argument (X : Complex_Matrix) return Real_Matrix
1007 renames Instantiations.Argument;
1010 (X : Complex_Matrix;
1011 Cycle : Real'Base) return Real_Matrix
1012 renames Instantiations.Argument;
1014 ----------------------------
1015 -- Compose_From_Cartesian --
1016 ----------------------------
1018 function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector
1019 renames Instantiations.Compose_From_Cartesian;
1021 function Compose_From_Cartesian
1023 Im : Real_Vector) return Complex_Vector
1024 renames Instantiations.Compose_From_Cartesian;
1026 function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix
1027 renames Instantiations.Compose_From_Cartesian;
1029 function Compose_From_Cartesian
1031 Im : Real_Matrix) return Complex_Matrix
1032 renames Instantiations.Compose_From_Cartesian;
1034 ------------------------
1035 -- Compose_From_Polar --
1036 ------------------------
1038 function Compose_From_Polar
1039 (Modulus : Real_Vector;
1040 Argument : Real_Vector) return Complex_Vector
1041 renames Instantiations.Compose_From_Polar;
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;
1049 function Compose_From_Polar
1050 (Modulus : Real_Matrix;
1051 Argument : Real_Matrix) return Complex_Matrix
1052 renames Instantiations.Compose_From_Polar;
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;
1064 function Conjugate (X : Complex_Vector) return Complex_Vector
1065 renames Instantiations.Conjugate;
1067 function Conjugate (X : Complex_Matrix) return Complex_Matrix
1068 renames Instantiations.Conjugate;
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;
1087 getrf (N, N, LU, N, Piv, Info'Access);
1090 raise Constraint_Error with "ill-conditioned matrix";
1094 Neg := Piv (1) /= 1;
1096 for J in 2 .. N loop
1097 Det := Det * LU (J, J);
1098 Neg := Neg xor (Piv (J) /= J);
1113 procedure Eigensystem
1114 (A : Complex_Matrix;
1115 Values : out Real_Vector;
1116 Vectors : out Complex_Matrix)
1118 Job_Z : aliased Character := 'V';
1119 Rng : aliased Character := 'A';
1120 Uplo : aliased Character := 'U';
1122 N : constant Natural := Length (A);
1123 W : BLAS_Real_Vector (Values'Range);
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;
1133 if Values'Length /= N then
1134 raise Constraint_Error with "wrong length for output vector";
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)
1142 raise Constraint_Error with "wrong dimensions for output matrix";
1149 -- Check for hermitian matrix ???
1150 -- Copy only required triangle ???
1154 -- Find size of work area
1157 (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1162 I_Supp_Z => I_Supp_Z,
1169 Info => Info'Access);
1172 raise Constraint_Error;
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));
1182 (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1187 I_Supp_Z => I_Supp_Z,
1189 L_Work => Work'Length,
1191 LR_Work => LR_Work'Length,
1193 LI_Work => LI_Work'Length,
1194 Info => Info'Access);
1197 raise Constraint_Error with "inverting non-Hermetian matrix";
1200 for J in Values'Range loop
1201 Values (J) := W (J);
1210 procedure Eigenvalues
1211 (A : Complex_Matrix;
1212 Values : out Real_Vector)
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);
1226 Info : aliased Integer;
1229 if Values'Length /= N then
1230 raise Constraint_Error with "wrong length for output vector";
1237 -- Check for hermitian matrix ???
1239 -- Find size of work area
1241 heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
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);
1253 raise Constraint_Error;
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));
1261 heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
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);
1273 raise Constraint_Error with "inverting singular matrix";
1276 for J in Values'Range loop
1277 Values (J) := W (J);
1282 function Eigenvalues (A : Complex_Matrix) return Real_Vector is
1283 R : Real_Vector (A'Range (1));
1293 function Im (X : Complex_Vector) return Real_Vector
1294 renames Instantiations.Im;
1296 function Im (X : Complex_Matrix) return Real_Matrix
1297 renames Instantiations.Im;
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;
1310 -- All computations are done using column-major order, but this works
1311 -- out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
1315 -- Compute LU decomposition
1322 Info => Info'Access);
1325 raise Constraint_Error with "inverting singular matrix";
1328 -- Determine size of work area
1336 Info => Info'Access);
1339 raise Constraint_Error;
1343 Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1346 -- Compute inverse from LU decomposition
1353 L_Work => Work'Length,
1354 Info => Info'Access);
1357 raise Constraint_Error with "inverting singular matrix";
1360 -- ??? Should iterate with gerfs, based on implementation advice
1364 function Inverse (A : Complex_Matrix) return Complex_Matrix is
1365 R : Complex_Matrix (A'Range (2), A'Range (1));
1375 function Modulus (X : Complex_Vector) return Real_Vector
1376 renames Instantiations.Modulus;
1378 function Modulus (X : Complex_Matrix) return Real_Matrix
1379 renames Instantiations.Modulus;
1385 function Re (X : Complex_Vector) return Real_Vector
1386 renames Instantiations.Re;
1388 function Re (X : Complex_Matrix) return Real_Matrix
1389 renames Instantiations.Re;
1396 (X : in out Complex_Matrix;
1398 renames Instantiations.Set_Im;
1401 (X : in out Complex_Vector;
1403 renames Instantiations.Set_Im;
1410 (X : in out Complex_Matrix;
1412 renames Instantiations.Set_Re;
1415 (X : in out Complex_Vector;
1417 renames Instantiations.Set_Re;
1424 (A : Complex_Matrix;
1426 B : out Complex_Vector)
1429 if Length (A) /= X'Length then
1430 raise Constraint_Error with
1431 "incompatible matrix and vector dimensions";
1434 -- ??? Should solve directly, is faster and more accurate
1436 B := Inverse (A) * X;
1440 (A : Complex_Matrix;
1442 B : out Complex_Matrix)
1445 if Length (A) /= X'Length (1) then
1446 raise Constraint_Error with "incompatible matrix dimensions";
1449 -- ??? Should solve directly, is faster and more accurate
1451 B := Inverse (A) * X;
1455 (A : Complex_Matrix;
1456 X : Complex_Vector) return Complex_Vector
1458 B : Complex_Vector (A'Range (2));
1464 function Solve (A, X : Complex_Matrix) return Complex_Matrix is
1465 B : Complex_Matrix (A'Range (2), X'Range (2));
1476 (X : Complex_Matrix) return Complex_Matrix
1478 R : Complex_Matrix (X'Range (2), X'Range (1));
1488 function Unit_Matrix
1490 First_1 : Integer := 1;
1491 First_2 : Integer := 1) return Complex_Matrix
1492 renames Instantiations.Unit_Matrix;
1498 function Unit_Vector
1501 First : Integer := 1) return Complex_Vector
1502 renames Instantiations.Unit_Vector;
1504 end Ada.Numerics.Generic_Complex_Arrays;