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-Hermitian 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
;