gcc/ada/a-ngcoar.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT COMPILER COMPONENTS                         --
   4 --                                                                          --
   5 --                   ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS                    --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --            Copyright (C) 2006-2009, Free Software Foundation, Inc.       --
  10 --                                                                          --
  11 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12 -- terms of the  GNU General Public License as published  by the Free Soft- --
  13 -- ware  Foundation;  either version 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.                                     --
  17 --                                                                          --
  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.               --
  21 --                                                                          --
  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/>.                                          --
  26 --                                                                          --
  27 -- GNAT was originally developed  by the GNAT team at  New York University. --
  28 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  29 --                                                                          --
  30 ------------------------------------------------------------------------------
  31
  32 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
  33 with System.Generic_Complex_BLAS;
  34 with System.Generic_Complex_LAPACK;
  35
  36 package body Ada.Numerics.Generic_Complex_Arrays is
  37
  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.
  42
  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.
  46
  47    --  Operations for solving linear systems and computing determinant,
  48    --  eigenvalues, eigensystem and inverse, are implemented using the
  49    --  LAPACK library.
  50
  51    type BLAS_Real_Vector is array (Integer range <>) of Real;
  52
  53    package BLAS is new System.Generic_Complex_BLAS
  54      (Real           => Real,
  55       Complex_Types  => Complex_Types,
  56       Complex_Vector => Complex_Vector,
  57       Complex_Matrix => Complex_Matrix);
  58
  59    package LAPACK is new System.Generic_Complex_LAPACK
  60      (Real           => Real,
  61       Real_Vector    => BLAS_Real_Vector,
  62       Complex_Types  => Complex_Types,
  63       Complex_Vector => Complex_Vector,
  64       Complex_Matrix => Complex_Matrix);
  65
  66    subtype Real is Real_Arrays.Real;
  67    --  Work around visibility bug ???
  68
  69    use BLAS, LAPACK;
  70
  71    --  Procedure versions of functions returning unconstrained values.
  72    --  This allows for inlining the function wrapper.
  73
  74    procedure Eigenvalues
  75      (A      : Complex_Matrix;
  76       Values : out Real_Vector);
  77
  78    procedure Inverse
  79      (A      : Complex_Matrix;
  80       R      : out Complex_Matrix);
  81
  82    procedure Solve
  83      (A      : Complex_Matrix;
  84       X      : Complex_Vector;
  85       B      : out Complex_Vector);
  86
  87    procedure Solve
  88      (A      : Complex_Matrix;
  89       X      : Complex_Matrix;
  90       B      : out Complex_Matrix);
  91
  92    procedure Transpose is new System.Generic_Array_Operations.Transpose
  93                                 (Scalar => Complex,
  94                                  Matrix => Complex_Matrix);
  95
  96    --  Helper function that raises a Constraint_Error is the argument is
  97    --  not a square matrix, and otherwise returns its length.
  98
  99    function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
 100
 101    --  Instantiating the following subprograms directly would lead to
 102    --  name clashes, so use a local package.
 103
 104    package Instantiations is
 105
 106       ---------
 107       -- "*" --
 108       ---------
 109
 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,
 116                              Operation     => "*");
 117
 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,
 124                              Operation     => "*");
 125
 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,
 132                              Operation     => "*");
 133
 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,
 140                              Operation     => "*");
 141
 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,
 148                              Zero          => (0.0, 0.0));
 149
 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,
 156                              Zero          => (0.0, 0.0));
 157
 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);
 165
 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);
 173
 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);
 181
 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,
 188                              Operation     => "*");
 189
 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,
 196                              Operation     => "*");
 197
 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,
 204                              Operation     => "*");
 205
 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,
 212                              Operation     => "*");
 213
 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,
 221                              Zero          => (0.0, 0.0));
 222
 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,
 230                              Zero          => (0.0, 0.0));
 231
 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,
 239                              Zero          => (0.0, 0.0));
 240
 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,
 248                              Zero          => (0.0, 0.0));
 249
 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,
 257                              Zero          => (0.0, 0.0));
 258
 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,
 266                              Zero          => (0.0, 0.0));
 267
 268       ---------
 269       -- "+" --
 270       ---------
 271
 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,
 277                              Operation     => "+");
 278
 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,
 286                              Operation     => "+");
 287
 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,
 295                              Operation     => "+");
 296
 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,
 304                              Operation     => "+");
 305
 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,
 311                              Operation     => "+");
 312
 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,
 320                              Operation     => "+");
 321
 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,
 329                              Operation     => "+");
 330
 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,
 338                              Operation     => "+");
 339
 340       ---------
 341       -- "-" --
 342       ---------
 343
 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,
 349                              Operation     => "-");
 350
 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,
 358                              Operation     => "-");
 359
 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,
 367                              Operation     => "-");
 368
 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,
 376                              Operation     => "-");
 377
 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,
 383                              Operation     => "-");
 384
 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,
 392                              Operation     => "-");
 393
 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,
 401                              Operation     => "-");
 402
 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,
 410                              Operation     => "-");
 411
 412       ---------
 413       -- "/" --
 414       ---------
 415
 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,
 422                              Operation     => "/");
 423
 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,
 430                              Operation     => "/");
 431
 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,
 438                              Operation     => "/");
 439
 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,
 446                              Operation     => "/");
 447
 448       --------------
 449       -- Argument --
 450       --------------
 451
 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);
 458
 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);
 466
 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);
 473
 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);
 481
 482       ----------------------------
 483       -- Compose_From_Cartesian --
 484       ----------------------------
 485
 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);
 492
 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);
 502
 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);
 509
 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);
 519
 520       ------------------------
 521       -- Compose_From_Polar --
 522       ------------------------
 523
 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);
 533
 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);
 544
 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);
 554
 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);
 565
 566       ---------------
 567       -- Conjugate --
 568       ---------------
 569
 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);
 576
 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);
 583
 584       --------
 585       -- Im --
 586       --------
 587
 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,
 593                              Operation     => Im);
 594
 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,
 600                              Operation     => Im);
 601
 602       -------------
 603       -- Modulus --
 604       -------------
 605
 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);
 612
 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);
 619
 620       --------
 621       -- Re --
 622       --------
 623
 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,
 629                              Operation     => Re);
 630
 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,
 636                              Operation     => Re);
 637
 638       ------------
 639       -- Set_Im --
 640       ------------
 641
 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,
 647                              Update        => Set_Im);
 648
 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,
 654                              Update        => Set_Im);
 655
 656       ------------
 657       -- Set_Re --
 658       ------------
 659
 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,
 665                              Update        => Set_Re);
 666
 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,
 672                              Update        => Set_Re);
 673
 674       -----------------
 675       -- Unit_Matrix --
 676       -----------------
 677
 678       function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix
 679                             (Scalar        => Complex,
 680                              Matrix        => Complex_Matrix,
 681                              Zero          => (0.0, 0.0),
 682                              One           => (1.0, 0.0));
 683
 684       function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector
 685                             (Scalar        => Complex,
 686                              Vector        => Complex_Vector,
 687                              Zero          => (0.0, 0.0),
 688                              One           => (1.0, 0.0));
 689
 690    end Instantiations;
 691
 692    ---------
 693    -- "*" --
 694    ---------
 695
 696    function "*"
 697      (Left  : Complex_Vector;
 698       Right : Complex_Vector) return Complex
 699    is
 700    begin
 701       if Left'Length /= Right'Length then
 702          raise Constraint_Error with
 703             "vectors are of different length in inner product";
 704       end if;
 705
 706       return dot (Left'Length, X => Left, Y => Right);
 707    end "*";
 708
 709    function "*"
 710      (Left  : Real_Vector;
 711       Right : Complex_Vector) return Complex
 712      renames Instantiations."*";
 713
 714    function "*"
 715      (Left  : Complex_Vector;
 716       Right : Real_Vector) return Complex
 717      renames Instantiations."*";
 718
 719    function "*"
 720      (Left  : Complex;
 721       Right : Complex_Vector) return Complex_Vector
 722      renames Instantiations."*";
 723
 724    function "*"
 725      (Left  : Complex_Vector;
 726       Right : Complex) return Complex_Vector
 727      renames Instantiations."*";
 728
 729    function "*"
 730      (Left  : Real'Base;
 731       Right : Complex_Vector) return Complex_Vector
 732      renames Instantiations."*";
 733
 734    function "*"
 735      (Left  : Complex_Vector;
 736       Right : Real'Base) return Complex_Vector
 737      renames Instantiations."*";
 738
 739    function "*"
 740      (Left  : Complex_Matrix;
 741       Right : Complex_Matrix)
 742       return  Complex_Matrix
 743    is
 744       R : Complex_Matrix (Left'Range (1), Right'Range (2));
 745
 746    begin
 747       if Left'Length (2) /= Right'Length (1) then
 748          raise Constraint_Error with
 749             "incompatible dimensions in matrix-matrix multiplication";
 750       end if;
 751
 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),
 757             A       => Right,
 758             Ld_A    => Right'Length (2),
 759             B       => Left,
 760             Ld_B    => Left'Length (2),
 761             C       => R,
 762             Ld_C    => R'Length (2));
 763
 764       return R;
 765    end "*";
 766
 767    function "*"
 768      (Left  : Complex_Vector;
 769       Right : Complex_Vector) return Complex_Matrix
 770      renames Instantiations."*";
 771
 772    function "*"
 773      (Left  : Complex_Vector;
 774       Right : Complex_Matrix) return Complex_Vector
 775    is
 776       R : Complex_Vector (Right'Range (2));
 777
 778    begin
 779       if Left'Length /= Right'Length (1) then
 780          raise Constraint_Error with
 781            "incompatible dimensions in vector-matrix multiplication";
 782       end if;
 783
 784       gemv (Trans => No_Trans'Access,
 785             M     => Right'Length (2),
 786             N     => Right'Length (1),
 787             A     => Right,
 788             Ld_A  => Right'Length (2),
 789             X     => Left,
 790             Y     => R);
 791
 792       return R;
 793    end "*";
 794
 795    function "*"
 796      (Left  : Complex_Matrix;
 797       Right : Complex_Vector) return Complex_Vector
 798    is
 799       R : Complex_Vector (Left'Range (1));
 800
 801    begin
 802       if Left'Length (2) /= Right'Length then
 803          raise Constraint_Error with
 804             "incompatible dimensions in matrix-vector multiplication";
 805       end if;
 806
 807       gemv (Trans => Trans'Access,
 808             M     => Left'Length (2),
 809             N     => Left'Length (1),
 810             A     => Left,
 811             Ld_A  => Left'Length (2),
 812             X     => Right,
 813             Y     => R);
 814
 815       return R;
 816    end "*";
 817
 818    function "*"
 819      (Left  : Real_Matrix;
 820       Right : Complex_Matrix) return Complex_Matrix
 821      renames Instantiations."*";
 822
 823    function "*"
 824      (Left  : Complex_Matrix;
 825       Right : Real_Matrix) return Complex_Matrix
 826      renames Instantiations."*";
 827
 828    function "*"
 829      (Left  : Real_Vector;
 830       Right : Complex_Vector) return Complex_Matrix
 831      renames Instantiations."*";
 832
 833    function "*"
 834      (Left  : Complex_Vector;
 835       Right : Real_Vector) return Complex_Matrix
 836      renames Instantiations."*";
 837
 838    function "*"
 839      (Left  : Real_Vector;
 840       Right : Complex_Matrix) return Complex_Vector
 841      renames Instantiations."*";
 842
 843    function "*"
 844      (Left  : Complex_Vector;
 845       Right : Real_Matrix) return Complex_Vector
 846      renames Instantiations."*";
 847
 848    function "*"
 849      (Left  : Real_Matrix;
 850       Right : Complex_Vector) return Complex_Vector
 851      renames Instantiations."*";
 852
 853    function "*"
 854      (Left  : Complex_Matrix;
 855       Right : Real_Vector) return Complex_Vector
 856      renames Instantiations."*";
 857
 858    function "*"
 859      (Left  : Complex;
 860       Right : Complex_Matrix) return Complex_Matrix
 861      renames Instantiations."*";
 862
 863    function "*"
 864      (Left  : Complex_Matrix;
 865       Right : Complex) return Complex_Matrix
 866      renames Instantiations."*";
 867
 868    function "*"
 869      (Left  : Real'Base;
 870       Right : Complex_Matrix) return Complex_Matrix
 871      renames Instantiations."*";
 872
 873    function "*"
 874      (Left  : Complex_Matrix;
 875       Right : Real'Base) return Complex_Matrix
 876      renames Instantiations."*";
 877
 878    ---------
 879    -- "+" --
 880    ---------
 881
 882    function "+" (Right : Complex_Vector) return Complex_Vector
 883      renames Instantiations."+";
 884
 885    function "+"
 886      (Left  : Complex_Vector;
 887       Right : Complex_Vector) return Complex_Vector
 888      renames Instantiations."+";
 889
 890    function "+"
 891      (Left  : Real_Vector;
 892       Right : Complex_Vector) return Complex_Vector
 893      renames Instantiations."+";
 894
 895    function "+"
 896      (Left  : Complex_Vector;
 897       Right : Real_Vector) return Complex_Vector
 898      renames Instantiations."+";
 899
 900    function "+" (Right : Complex_Matrix) return Complex_Matrix
 901      renames Instantiations."+";
 902
 903    function "+"
 904      (Left  : Complex_Matrix;
 905       Right : Complex_Matrix) return Complex_Matrix
 906      renames Instantiations."+";
 907
 908    function "+"
 909      (Left  : Real_Matrix;
 910       Right : Complex_Matrix) return Complex_Matrix
 911      renames Instantiations."+";
 912
 913    function "+"
 914      (Left  : Complex_Matrix;
 915       Right : Real_Matrix) return Complex_Matrix
 916      renames Instantiations."+";
 917
 918    ---------
 919    -- "-" --
 920    ---------
 921
 922    function "-"
 923      (Right : Complex_Vector) return Complex_Vector
 924      renames Instantiations."-";
 925
 926    function "-"
 927      (Left  : Complex_Vector;
 928       Right : Complex_Vector) return Complex_Vector
 929      renames Instantiations."-";
 930
 931    function "-"
 932      (Left  : Real_Vector;
 933       Right : Complex_Vector) return Complex_Vector
 934       renames Instantiations."-";
 935
 936    function "-"
 937      (Left  : Complex_Vector;
 938       Right : Real_Vector) return Complex_Vector
 939      renames Instantiations."-";
 940
 941    function "-" (Right : Complex_Matrix) return Complex_Matrix
 942      renames Instantiations."-";
 943
 944    function "-"
 945      (Left  : Complex_Matrix;
 946       Right : Complex_Matrix) return Complex_Matrix
 947      renames Instantiations."-";
 948
 949    function "-"
 950      (Left  : Real_Matrix;
 951       Right : Complex_Matrix) return Complex_Matrix
 952      renames Instantiations."-";
 953
 954    function "-"
 955      (Left  : Complex_Matrix;
 956       Right : Real_Matrix) return Complex_Matrix
 957      renames Instantiations."-";
 958
 959    ---------
 960    -- "/" --
 961    ---------
 962
 963    function "/"
 964      (Left  : Complex_Vector;
 965       Right : Complex) return Complex_Vector
 966      renames Instantiations."/";
 967
 968    function "/"
 969      (Left  : Complex_Vector;
 970       Right : Real'Base) return Complex_Vector
 971      renames Instantiations."/";
 972
 973    function "/"
 974      (Left  : Complex_Matrix;
 975       Right : Complex) return Complex_Matrix
 976      renames Instantiations."/";
 977
 978    function "/"
 979      (Left  : Complex_Matrix;
 980       Right : Real'Base) return Complex_Matrix
 981      renames Instantiations."/";
 982
 983    -----------
 984    -- "abs" --
 985    -----------
 986
 987    function "abs" (Right : Complex_Vector) return Complex is
 988    begin
 989       return (nrm2 (Right'Length, Right), 0.0);
 990    end "abs";
 991
 992    --------------
 993    -- Argument --
 994    --------------
 995
 996    function Argument (X : Complex_Vector) return Real_Vector
 997      renames Instantiations.Argument;
 998
 999    function Argument
1000      (X     : Complex_Vector;
1001       Cycle : Real'Base) return Real_Vector
1002      renames Instantiations.Argument;
1003
1004    function Argument (X : Complex_Matrix) return Real_Matrix
1005      renames Instantiations.Argument;
1006
1007    function Argument
1008      (X     : Complex_Matrix;
1009       Cycle : Real'Base) return Real_Matrix
1010      renames Instantiations.Argument;
1011
1012    ----------------------------
1013    -- Compose_From_Cartesian --
1014    ----------------------------
1015
1016    function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector
1017      renames Instantiations.Compose_From_Cartesian;
1018
1019    function Compose_From_Cartesian
1020      (Re : Real_Vector;
1021       Im : Real_Vector) return Complex_Vector
1022      renames Instantiations.Compose_From_Cartesian;
1023
1024    function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix
1025      renames Instantiations.Compose_From_Cartesian;
1026
1027    function Compose_From_Cartesian
1028      (Re : Real_Matrix;
1029       Im : Real_Matrix) return Complex_Matrix
1030      renames Instantiations.Compose_From_Cartesian;
1031
1032    ------------------------
1033    -- Compose_From_Polar --
1034    ------------------------
1035
1036    function Compose_From_Polar
1037      (Modulus  : Real_Vector;
1038       Argument : Real_Vector) return Complex_Vector
1039      renames Instantiations.Compose_From_Polar;
1040
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;
1046
1047    function Compose_From_Polar
1048      (Modulus  : Real_Matrix;
1049       Argument : Real_Matrix) return Complex_Matrix
1050      renames Instantiations.Compose_From_Polar;
1051
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;
1057
1058    ---------------
1059    -- Conjugate --
1060    ---------------
1061
1062    function Conjugate (X : Complex_Vector) return Complex_Vector
1063      renames Instantiations.Conjugate;
1064
1065    function Conjugate (X : Complex_Matrix) return Complex_Matrix
1066      renames Instantiations.Conjugate;
1067
1068    -----------------
1069    -- Determinant --
1070    -----------------
1071
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;
1077       Neg  : Boolean;
1078       Det  : Complex;
1079
1080    begin
1081       if N = 0 then
1082          return (0.0, 0.0);
1083       end if;
1084
1085       getrf (N, N, LU, N, Piv, Info'Access);
1086
1087       if Info /= 0 then
1088          raise Constraint_Error with "ill-conditioned matrix";
1089       end if;
1090
1091       Det := LU (1, 1);
1092       Neg := Piv (1) /= 1;
1093
1094       for J in 2 .. N loop
1095          Det := Det * LU (J, J);
1096          Neg := Neg xor (Piv (J) /= J);
1097       end loop;
1098
1099       if Neg then
1100          return -Det;
1101
1102       else
1103          return Det;
1104       end if;
1105    end Determinant;
1106
1107    -----------------
1108    -- Eigensystem --
1109    -----------------
1110
1111    procedure Eigensystem
1112      (A       : Complex_Matrix;
1113       Values  : out Real_Vector;
1114       Vectors : out Complex_Matrix)
1115    is
1116       Job_Z    : aliased Character := 'V';
1117       Rng      : aliased Character := 'A';
1118       Uplo     : aliased Character := 'U';
1119
1120       N        : constant Natural := Length (A);
1121       W        : BLAS_Real_Vector (Values'Range);
1122       M        : Integer;
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;
1129
1130    begin
1131       if Values'Length /= N then
1132          raise Constraint_Error with "wrong length for output vector";
1133       end if;
1134
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)
1139       then
1140          raise Constraint_Error with "wrong dimensions for output matrix";
1141       end if;
1142
1143       if N = 0 then
1144          return;
1145       end if;
1146
1147       --  Check for hermitian matrix ???
1148       --  Copy only required triangle ???
1149
1150       B := A;
1151
1152       --  Find size of work area
1153
1154       heevr
1155         (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1156          M        => M,
1157          W        => W,
1158          Z        => Vectors,
1159          Ld_Z     => N,
1160          I_Supp_Z => I_Supp_Z,
1161          Work     => L_Work,
1162          L_Work   => -1,
1163          R_Work   => LR_Work,
1164          LR_Work  => -1,
1165          I_Work   => LI_Work,
1166          LI_Work  => -1,
1167          Info     => Info'Access);
1168
1169       if Info /= 0 then
1170          raise Constraint_Error;
1171       end if;
1172
1173       declare
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));
1177
1178       begin
1179          heevr
1180            (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1181             M        => M,
1182             W        => W,
1183             Z        => Vectors,
1184             Ld_Z     => N,
1185             I_Supp_Z => I_Supp_Z,
1186             Work     => Work,
1187             L_Work   => Work'Length,
1188             R_Work   => R_Work,
1189             LR_Work  => LR_Work'Length,
1190             I_Work   => I_Work,
1191             LI_Work  => LI_Work'Length,
1192             Info     => Info'Access);
1193
1194          if Info /= 0 then
1195             raise Constraint_Error with "inverting non-Hermitian matrix";
1196          end if;
1197
1198          for J in Values'Range loop
1199             Values (J) := W (J);
1200          end loop;
1201       end;
1202    end Eigensystem;
1203
1204    -----------------
1205    -- Eigenvalues --
1206    -----------------
1207
1208    procedure Eigenvalues
1209      (A      : Complex_Matrix;
1210       Values : out Real_Vector)
1211    is
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);
1223       M        : Integer;
1224       Info     : aliased Integer;
1225
1226    begin
1227       if Values'Length /= N then
1228          raise Constraint_Error with "wrong length for output vector";
1229       end if;
1230
1231       if N = 0 then
1232          return;
1233       end if;
1234
1235       --  Check for hermitian matrix ???
1236
1237       --  Find size of work area
1238
1239       heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1240              M        => M,
1241              W        => W,
1242              Z        => Z,
1243              Ld_Z     => 1,
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);
1249
1250       if Info /= 0 then
1251          raise Constraint_Error;
1252       end if;
1253
1254       declare
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));
1258       begin
1259          heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1260                 M        => M,
1261                 W        => W,
1262                 Z        => Z,
1263                 Ld_Z     => 1,
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);
1269
1270          if Info /= 0 then
1271             raise Constraint_Error with "inverting singular matrix";
1272          end if;
1273
1274          for J in Values'Range loop
1275             Values (J) := W (J);
1276          end loop;
1277       end;
1278    end Eigenvalues;
1279
1280    function Eigenvalues (A : Complex_Matrix) return Real_Vector is
1281       R : Real_Vector (A'Range (1));
1282    begin
1283       Eigenvalues (A, R);
1284       return R;
1285    end Eigenvalues;
1286
1287    --------
1288    -- Im --
1289    --------
1290
1291    function Im (X : Complex_Vector) return Real_Vector
1292      renames Instantiations.Im;
1293
1294    function Im (X : Complex_Matrix) return Real_Matrix
1295      renames Instantiations.Im;
1296
1297    -------------
1298    -- Inverse --
1299    -------------
1300
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;
1306
1307    begin
1308       --  All computations are done using column-major order, but this works
1309       --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
1310
1311       R := A;
1312
1313       --  Compute LU decomposition
1314
1315       getrf (M      => N,
1316              N      => N,
1317              A      => R,
1318              Ld_A   => N,
1319              I_Piv  => Piv,
1320              Info   => Info'Access);
1321
1322       if Info /= 0 then
1323          raise Constraint_Error with "inverting singular matrix";
1324       end if;
1325
1326       --  Determine size of work area
1327
1328       getri (N      => N,
1329              A      => R,
1330              Ld_A   => N,
1331              I_Piv  => Piv,
1332              Work   => L_Work,
1333              L_Work => -1,
1334              Info   => Info'Access);
1335
1336       if Info /= 0 then
1337          raise Constraint_Error;
1338       end if;
1339
1340       declare
1341          Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1342
1343       begin
1344          --  Compute inverse from LU decomposition
1345
1346          getri (N      => N,
1347                 A      => R,
1348                 Ld_A   => N,
1349                 I_Piv  => Piv,
1350                 Work   => Work,
1351                 L_Work => Work'Length,
1352                 Info   => Info'Access);
1353
1354          if Info /= 0 then
1355             raise Constraint_Error with "inverting singular matrix";
1356          end if;
1357
1358          --  ??? Should iterate with gerfs, based on implementation advice
1359       end;
1360    end Inverse;
1361
1362    function Inverse (A : Complex_Matrix) return Complex_Matrix is
1363       R : Complex_Matrix (A'Range (2), A'Range (1));
1364    begin
1365       Inverse (A, R);
1366       return R;
1367    end Inverse;
1368
1369    -------------
1370    -- Modulus --
1371    -------------
1372
1373    function Modulus (X : Complex_Vector) return Real_Vector
1374      renames Instantiations.Modulus;
1375
1376    function Modulus (X : Complex_Matrix) return Real_Matrix
1377      renames Instantiations.Modulus;
1378
1379    --------
1380    -- Re --
1381    --------
1382
1383    function Re (X : Complex_Vector) return Real_Vector
1384      renames Instantiations.Re;
1385
1386    function Re (X : Complex_Matrix) return Real_Matrix
1387      renames Instantiations.Re;
1388
1389    ------------
1390    -- Set_Im --
1391    ------------
1392
1393    procedure Set_Im
1394      (X  : in out Complex_Matrix;
1395       Im : Real_Matrix)
1396      renames Instantiations.Set_Im;
1397
1398    procedure Set_Im
1399      (X  : in out Complex_Vector;
1400       Im : Real_Vector)
1401      renames Instantiations.Set_Im;
1402
1403    ------------
1404    -- Set_Re --
1405    ------------
1406
1407    procedure Set_Re
1408      (X  : in out Complex_Matrix;
1409       Re : Real_Matrix)
1410      renames Instantiations.Set_Re;
1411
1412    procedure Set_Re
1413      (X  : in out Complex_Vector;
1414       Re : Real_Vector)
1415      renames Instantiations.Set_Re;
1416
1417    -----------
1418    -- Solve --
1419    -----------
1420
1421    procedure Solve
1422      (A : Complex_Matrix;
1423       X : Complex_Vector;
1424       B : out Complex_Vector)
1425    is
1426    begin
1427       if Length (A) /= X'Length then
1428          raise Constraint_Error with
1429            "incompatible matrix and vector dimensions";
1430       end if;
1431
1432       --  ??? Should solve directly, is faster and more accurate
1433
1434       B := Inverse (A) * X;
1435    end Solve;
1436
1437    procedure Solve
1438      (A : Complex_Matrix;
1439       X : Complex_Matrix;
1440       B : out Complex_Matrix)
1441    is
1442    begin
1443       if Length (A) /= X'Length (1) then
1444          raise Constraint_Error with "incompatible matrix dimensions";
1445       end if;
1446
1447       --  ??? Should solve directly, is faster and more accurate
1448
1449       B := Inverse (A) * X;
1450    end Solve;
1451
1452    function Solve
1453      (A : Complex_Matrix;
1454       X : Complex_Vector) return Complex_Vector
1455    is
1456       B : Complex_Vector (A'Range (2));
1457    begin
1458       Solve (A, X, B);
1459       return B;
1460    end Solve;
1461
1462    function Solve (A, X : Complex_Matrix) return Complex_Matrix is
1463       B : Complex_Matrix (A'Range (2), X'Range (2));
1464    begin
1465       Solve (A, X, B);
1466       return B;
1467    end Solve;
1468
1469    ---------------
1470    -- Transpose --
1471    ---------------
1472
1473    function Transpose
1474      (X : Complex_Matrix) return Complex_Matrix
1475    is
1476       R : Complex_Matrix (X'Range (2), X'Range (1));
1477    begin
1478       Transpose (X, R);
1479       return R;
1480    end Transpose;
1481
1482    -----------------
1483    -- Unit_Matrix --
1484    -----------------
1485
1486    function Unit_Matrix
1487      (Order   : Positive;
1488       First_1 : Integer := 1;
1489       First_2 : Integer := 1) return Complex_Matrix
1490      renames Instantiations.Unit_Matrix;
1491
1492    -----------------
1493    -- Unit_Vector --
1494    -----------------
1495
1496    function Unit_Vector
1497      (Index : Integer;
1498       Order : Positive;
1499       First : Integer := 1) return Complex_Vector
1500      renames Instantiations.Unit_Vector;
1501
1502 end Ada.Numerics.Generic_Complex_Arrays;