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
;