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, 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 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.                                              --
  21 --                                                                          --
  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.                                      --
  28 --                                                                          --
  29 -- GNAT was originally developed  by the GNAT team at  New York University. --
  30 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  31 --                                                                          --
  32 ------------------------------------------------------------------------------
  33
  34 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
  35 with System.Generic_Complex_BLAS;
  36 with System.Generic_Complex_LAPACK;
  37
  38 package body Ada.Numerics.Generic_Complex_Arrays is
  39
  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.
  44
  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.
  48
  49    --  Operations for solving linear systems and computing determinant,
  50    --  eigenvalues, eigensystem and inverse, are implemented using the
  51    --  LAPACK library.
  52
  53    type BLAS_Real_Vector is array (Integer range <>) of Real;
  54
  55    package BLAS is new System.Generic_Complex_BLAS
  56      (Real           => Real,
  57       Complex_Types  => Complex_Types,
  58       Complex_Vector => Complex_Vector,
  59       Complex_Matrix => Complex_Matrix);
  60
  61    package LAPACK is new System.Generic_Complex_LAPACK
  62      (Real           => Real,
  63       Real_Vector    => BLAS_Real_Vector,
  64       Complex_Types  => Complex_Types,
  65       Complex_Vector => Complex_Vector,
  66       Complex_Matrix => Complex_Matrix);
  67
  68    use BLAS, LAPACK;
  69
  70    --  Procedure versions of functions returning unconstrained values.
  71    --  This allows for inlining the function wrapper.
  72
  73    procedure Eigenvalues
  74      (A      : Complex_Matrix;
  75       Values : out Real_Vector);
  76
  77    procedure Inverse
  78      (A      : Complex_Matrix;
  79       R      : out Complex_Matrix);
  80
  81    procedure Solve
  82      (A      : Complex_Matrix;
  83       X      : Complex_Vector;
  84       B      : out Complex_Vector);
  85
  86    procedure Solve
  87      (A      : Complex_Matrix;
  88       X      : Complex_Matrix;
  89       B      : out Complex_Matrix);
  90
  91    procedure Transpose is new System.Generic_Array_Operations.Transpose
  92                                 (Scalar => Complex,
  93                                  Matrix => Complex_Matrix);
  94
  95    --  Helper function that raises a Constraint_Error is the argument is
  96    --  not a square matrix, and otherwise returns its length.
  97
  98    function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
  99
 100    --  Instantiating the following subprograms directly would lead to
 101    --  name clashes, so use a local package.
 102
 103    package Instantiations is
 104
 105       ---------
 106       -- "*" --
 107       ---------
 108
 109       function "*" is new Vector_Scalar_Elementwise_Operation
 110                             (Left_Scalar   => Complex,
 111                              Right_Scalar  => Complex,
 112                              Result_Scalar => Complex,
 113                              Left_Vector   => Complex_Vector,
 114                              Result_Vector => Complex_Vector,
 115                              Operation     => "*");
 116
 117       function "*" is new Vector_Scalar_Elementwise_Operation
 118                             (Left_Scalar   => Complex,
 119                              Right_Scalar  => Real'Base,
 120                              Result_Scalar => Complex,
 121                              Left_Vector   => Complex_Vector,
 122                              Result_Vector => Complex_Vector,
 123                              Operation     => "*");
 124
 125       function "*" is new Scalar_Vector_Elementwise_Operation
 126                             (Left_Scalar   => Complex,
 127                              Right_Scalar  => Complex,
 128                              Result_Scalar => Complex,
 129                              Right_Vector  => Complex_Vector,
 130                              Result_Vector => Complex_Vector,
 131                              Operation     => "*");
 132
 133       function "*" is new Scalar_Vector_Elementwise_Operation
 134                             (Left_Scalar   => Real'Base,
 135                              Right_Scalar  => Complex,
 136                              Result_Scalar => Complex,
 137                              Right_Vector  => Complex_Vector,
 138                              Result_Vector => Complex_Vector,
 139                              Operation     => "*");
 140
 141       function "*" is new Inner_Product
 142                             (Left_Scalar   => Complex,
 143                              Right_Scalar  => Real'Base,
 144                              Result_Scalar => Complex,
 145                              Left_Vector   => Complex_Vector,
 146                              Right_Vector  => Real_Vector,
 147                              Zero          => (0.0, 0.0));
 148
 149       function "*" is new Inner_Product
 150                             (Left_Scalar   => Real'Base,
 151                              Right_Scalar  => Complex,
 152                              Result_Scalar => Complex,
 153                              Left_Vector   => Real_Vector,
 154                              Right_Vector  => Complex_Vector,
 155                              Zero          => (0.0, 0.0));
 156
 157       function "*" is new Outer_Product
 158                             (Left_Scalar   => Complex,
 159                              Right_Scalar  => Complex,
 160                              Result_Scalar => Complex,
 161                              Left_Vector   => Complex_Vector,
 162                              Right_Vector  => Complex_Vector,
 163                              Matrix        => Complex_Matrix);
 164
 165       function "*" is new Outer_Product
 166                             (Left_Scalar   => Real'Base,
 167                              Right_Scalar  => Complex,
 168                              Result_Scalar => Complex,
 169                              Left_Vector   => Real_Vector,
 170                              Right_Vector  => Complex_Vector,
 171                              Matrix        => Complex_Matrix);
 172
 173       function "*" is new Outer_Product
 174                             (Left_Scalar   => Complex,
 175                              Right_Scalar  => Real'Base,
 176                              Result_Scalar => Complex,
 177                              Left_Vector   => Complex_Vector,
 178                              Right_Vector  => Real_Vector,
 179                              Matrix        => Complex_Matrix);
 180
 181       function "*" is new Matrix_Scalar_Elementwise_Operation
 182                             (Left_Scalar   => Complex,
 183                              Right_Scalar  => Complex,
 184                              Result_Scalar => Complex,
 185                              Left_Matrix   => Complex_Matrix,
 186                              Result_Matrix => Complex_Matrix,
 187                              Operation     => "*");
 188
 189       function "*" is new Matrix_Scalar_Elementwise_Operation
 190                             (Left_Scalar   => Complex,
 191                              Right_Scalar  => Real'Base,
 192                              Result_Scalar => Complex,
 193                              Left_Matrix   => Complex_Matrix,
 194                              Result_Matrix => Complex_Matrix,
 195                              Operation     => "*");
 196
 197       function "*" is new Scalar_Matrix_Elementwise_Operation
 198                             (Left_Scalar   => Complex,
 199                              Right_Scalar  => Complex,
 200                              Result_Scalar => Complex,
 201                              Right_Matrix  => Complex_Matrix,
 202                              Result_Matrix => Complex_Matrix,
 203                              Operation     => "*");
 204
 205       function "*" is new Scalar_Matrix_Elementwise_Operation
 206                             (Left_Scalar   => Real'Base,
 207                              Right_Scalar  => Complex,
 208                              Result_Scalar => Complex,
 209                              Right_Matrix  => Complex_Matrix,
 210                              Result_Matrix => Complex_Matrix,
 211                              Operation     => "*");
 212
 213       function "*" is new Matrix_Vector_Product
 214                             (Left_Scalar   => Real'Base,
 215                              Right_Scalar  => Complex,
 216                              Result_Scalar => Complex,
 217                              Matrix        => Real_Matrix,
 218                              Right_Vector  => Complex_Vector,
 219                              Result_Vector => Complex_Vector,
 220                              Zero          => (0.0, 0.0));
 221
 222       function "*" is new Matrix_Vector_Product
 223                             (Left_Scalar   => Complex,
 224                              Right_Scalar  => Real'Base,
 225                              Result_Scalar => Complex,
 226                              Matrix        => Complex_Matrix,
 227                              Right_Vector  => Real_Vector,
 228                              Result_Vector => Complex_Vector,
 229                              Zero          => (0.0, 0.0));
 230
 231       function "*" is new Vector_Matrix_Product
 232                             (Left_Scalar   => Real'Base,
 233                              Right_Scalar  => Complex,
 234                              Result_Scalar => Complex,
 235                              Left_Vector   => Real_Vector,
 236                              Matrix        => Complex_Matrix,
 237                              Result_Vector => Complex_Vector,
 238                              Zero          => (0.0, 0.0));
 239
 240       function "*" is new Vector_Matrix_Product
 241                             (Left_Scalar   => Complex,
 242                              Right_Scalar  => Real'Base,
 243                              Result_Scalar => Complex,
 244                              Left_Vector   => Complex_Vector,
 245                              Matrix        => Real_Matrix,
 246                              Result_Vector => Complex_Vector,
 247                              Zero          => (0.0, 0.0));
 248
 249       function "*" is new Matrix_Matrix_Product
 250                             (Left_Scalar   => Real'Base,
 251                              Right_Scalar  => Complex,
 252                              Result_Scalar => Complex,
 253                              Left_Matrix   => Real_Matrix,
 254                              Right_Matrix  => Complex_Matrix,
 255                              Result_Matrix => Complex_Matrix,
 256                              Zero          => (0.0, 0.0));
 257
 258       function "*" is new Matrix_Matrix_Product
 259                             (Left_Scalar   => Complex,
 260                              Right_Scalar  => Real'Base,
 261                              Result_Scalar => Complex,
 262                              Left_Matrix   => Complex_Matrix,
 263                              Right_Matrix  => Real_Matrix,
 264                              Result_Matrix => Complex_Matrix,
 265                              Zero          => (0.0, 0.0));
 266
 267       ---------
 268       -- "+" --
 269       ---------
 270
 271       function "+" is new Vector_Elementwise_Operation
 272                             (X_Scalar      => Complex,
 273                              Result_Scalar => Complex,
 274                              X_Vector      => Complex_Vector,
 275                              Result_Vector => Complex_Vector,
 276                              Operation     => "+");
 277
 278       function "+" is new Vector_Vector_Elementwise_Operation
 279                             (Left_Scalar   => Complex,
 280                              Right_Scalar  => Complex,
 281                              Result_Scalar => Complex,
 282                              Left_Vector   => Complex_Vector,
 283                              Right_Vector  => Complex_Vector,
 284                              Result_Vector => Complex_Vector,
 285                              Operation     => "+");
 286
 287       function "+" is new Vector_Vector_Elementwise_Operation
 288                             (Left_Scalar   => Real'Base,
 289                              Right_Scalar  => Complex,
 290                              Result_Scalar => Complex,
 291                              Left_Vector   => Real_Vector,
 292                              Right_Vector  => Complex_Vector,
 293                              Result_Vector => Complex_Vector,
 294                              Operation     => "+");
 295
 296       function "+" is new Vector_Vector_Elementwise_Operation
 297                             (Left_Scalar   => Complex,
 298                              Right_Scalar  => Real'Base,
 299                              Result_Scalar => Complex,
 300                              Left_Vector   => Complex_Vector,
 301                              Right_Vector  => Real_Vector,
 302                              Result_Vector => Complex_Vector,
 303                              Operation     => "+");
 304
 305       function "+" is new Matrix_Elementwise_Operation
 306                             (X_Scalar      => Complex,
 307                              Result_Scalar => Complex,
 308                              X_Matrix      => Complex_Matrix,
 309                              Result_Matrix => Complex_Matrix,
 310                              Operation     => "+");
 311
 312       function "+" is new Matrix_Matrix_Elementwise_Operation
 313                             (Left_Scalar   => Complex,
 314                              Right_Scalar  => Complex,
 315                              Result_Scalar => Complex,
 316                              Left_Matrix   => Complex_Matrix,
 317                              Right_Matrix  => Complex_Matrix,
 318                              Result_Matrix => Complex_Matrix,
 319                              Operation     => "+");
 320
 321       function "+" is new Matrix_Matrix_Elementwise_Operation
 322                             (Left_Scalar   => Real'Base,
 323                              Right_Scalar  => Complex,
 324                              Result_Scalar => Complex,
 325                              Left_Matrix   => Real_Matrix,
 326                              Right_Matrix  => Complex_Matrix,
 327                              Result_Matrix => Complex_Matrix,
 328                              Operation     => "+");
 329
 330       function "+" is new Matrix_Matrix_Elementwise_Operation
 331                             (Left_Scalar   => Complex,
 332                              Right_Scalar  => Real'Base,
 333                              Result_Scalar => Complex,
 334                              Left_Matrix   => Complex_Matrix,
 335                              Right_Matrix  => Real_Matrix,
 336                              Result_Matrix => Complex_Matrix,
 337                              Operation     => "+");
 338
 339       ---------
 340       -- "-" --
 341       ---------
 342
 343       function "-" is new Vector_Elementwise_Operation
 344                             (X_Scalar      => Complex,
 345                              Result_Scalar => Complex,
 346                              X_Vector      => Complex_Vector,
 347                              Result_Vector => Complex_Vector,
 348                              Operation     => "-");
 349
 350       function "-" is new Vector_Vector_Elementwise_Operation
 351                             (Left_Scalar   => Complex,
 352                              Right_Scalar  => Complex,
 353                              Result_Scalar => Complex,
 354                              Left_Vector   => Complex_Vector,
 355                              Right_Vector  => Complex_Vector,
 356                              Result_Vector => Complex_Vector,
 357                              Operation     => "-");
 358
 359       function "-" is new Vector_Vector_Elementwise_Operation
 360                             (Left_Scalar   => Real'Base,
 361                              Right_Scalar  => Complex,
 362                              Result_Scalar => Complex,
 363                              Left_Vector   => Real_Vector,
 364                              Right_Vector  => Complex_Vector,
 365                              Result_Vector => Complex_Vector,
 366                              Operation     => "-");
 367
 368       function "-" is new Vector_Vector_Elementwise_Operation
 369                             (Left_Scalar   => Complex,
 370                              Right_Scalar  => Real'Base,
 371                              Result_Scalar => Complex,
 372                              Left_Vector   => Complex_Vector,
 373                              Right_Vector  => Real_Vector,
 374                              Result_Vector => Complex_Vector,
 375                              Operation     => "-");
 376
 377       function "-" is new Matrix_Elementwise_Operation
 378                             (X_Scalar      => Complex,
 379                              Result_Scalar => Complex,
 380                              X_Matrix      => Complex_Matrix,
 381                              Result_Matrix => Complex_Matrix,
 382                              Operation     => "-");
 383
 384       function "-" is new Matrix_Matrix_Elementwise_Operation
 385                             (Left_Scalar   => Complex,
 386                              Right_Scalar  => Complex,
 387                              Result_Scalar => Complex,
 388                              Left_Matrix   => Complex_Matrix,
 389                              Right_Matrix  => Complex_Matrix,
 390                              Result_Matrix => Complex_Matrix,
 391                              Operation     => "-");
 392
 393       function "-" is new Matrix_Matrix_Elementwise_Operation
 394                             (Left_Scalar   => Real'Base,
 395                              Right_Scalar  => Complex,
 396                              Result_Scalar => Complex,
 397                              Left_Matrix   => Real_Matrix,
 398                              Right_Matrix  => Complex_Matrix,
 399                              Result_Matrix => Complex_Matrix,
 400                              Operation     => "-");
 401
 402       function "-" is new Matrix_Matrix_Elementwise_Operation
 403                             (Left_Scalar   => Complex,
 404                              Right_Scalar  => Real'Base,
 405                              Result_Scalar => Complex,
 406                              Left_Matrix   => Complex_Matrix,
 407                              Right_Matrix  => Real_Matrix,
 408                              Result_Matrix => Complex_Matrix,
 409                              Operation     => "-");
 410
 411       ---------
 412       -- "/" --
 413       ---------
 414
 415       function "/" is new Vector_Scalar_Elementwise_Operation
 416                             (Left_Scalar   => Complex,
 417                              Right_Scalar  => Complex,
 418                              Result_Scalar => Complex,
 419                              Left_Vector   => Complex_Vector,
 420                              Result_Vector => Complex_Vector,
 421                              Operation     => "/");
 422
 423       function "/" is new Vector_Scalar_Elementwise_Operation
 424                             (Left_Scalar   => Complex,
 425                              Right_Scalar  => Real'Base,
 426                              Result_Scalar => Complex,
 427                              Left_Vector   => Complex_Vector,
 428                              Result_Vector => Complex_Vector,
 429                              Operation     => "/");
 430
 431       function "/" is new Matrix_Scalar_Elementwise_Operation
 432                             (Left_Scalar   => Complex,
 433                              Right_Scalar  => Complex,
 434                              Result_Scalar => Complex,
 435                              Left_Matrix   => Complex_Matrix,
 436                              Result_Matrix => Complex_Matrix,
 437                              Operation     => "/");
 438
 439       function "/" is new Matrix_Scalar_Elementwise_Operation
 440                             (Left_Scalar   => Complex,
 441                              Right_Scalar  => Real'Base,
 442                              Result_Scalar => Complex,
 443                              Left_Matrix   => Complex_Matrix,
 444                              Result_Matrix => Complex_Matrix,
 445                              Operation     => "/");
 446
 447       --------------
 448       -- Argument --
 449       --------------
 450
 451       function Argument is new Vector_Elementwise_Operation
 452                             (X_Scalar      => Complex,
 453                              Result_Scalar => Real'Base,
 454                              X_Vector      => Complex_Vector,
 455                              Result_Vector => Real_Vector,
 456                              Operation     => Argument);
 457
 458       function Argument is new Vector_Scalar_Elementwise_Operation
 459                             (Left_Scalar   => Complex,
 460                              Right_Scalar  => Real'Base,
 461                              Result_Scalar => Real'Base,
 462                              Left_Vector   => Complex_Vector,
 463                              Result_Vector => Real_Vector,
 464                              Operation     => Argument);
 465
 466       function Argument is new Matrix_Elementwise_Operation
 467                             (X_Scalar      => Complex,
 468                              Result_Scalar => Real'Base,
 469                              X_Matrix      => Complex_Matrix,
 470                              Result_Matrix => Real_Matrix,
 471                              Operation     => Argument);
 472
 473       function Argument is new Matrix_Scalar_Elementwise_Operation
 474                             (Left_Scalar   => Complex,
 475                              Right_Scalar  => Real'Base,
 476                              Result_Scalar => Real'Base,
 477                              Left_Matrix   => Complex_Matrix,
 478                              Result_Matrix => Real_Matrix,
 479                              Operation     => Argument);
 480
 481       ----------------------------
 482       -- Compose_From_Cartesian --
 483       ----------------------------
 484
 485       function Compose_From_Cartesian is new Vector_Elementwise_Operation
 486                             (X_Scalar      => Real'Base,
 487                              Result_Scalar => Complex,
 488                              X_Vector      => Real_Vector,
 489                              Result_Vector => Complex_Vector,
 490                              Operation     => Compose_From_Cartesian);
 491
 492       function Compose_From_Cartesian is
 493          new Vector_Vector_Elementwise_Operation
 494                             (Left_Scalar   => Real'Base,
 495                              Right_Scalar  => Real'Base,
 496                              Result_Scalar => Complex,
 497                              Left_Vector   => Real_Vector,
 498                              Right_Vector  => Real_Vector,
 499                              Result_Vector => Complex_Vector,
 500                              Operation     => Compose_From_Cartesian);
 501
 502       function Compose_From_Cartesian is new Matrix_Elementwise_Operation
 503                             (X_Scalar      => Real'Base,
 504                              Result_Scalar => Complex,
 505                              X_Matrix      => Real_Matrix,
 506                              Result_Matrix => Complex_Matrix,
 507                              Operation     => Compose_From_Cartesian);
 508
 509       function Compose_From_Cartesian is
 510          new Matrix_Matrix_Elementwise_Operation
 511                             (Left_Scalar   => Real'Base,
 512                              Right_Scalar  => Real'Base,
 513                              Result_Scalar => Complex,
 514                              Left_Matrix   => Real_Matrix,
 515                              Right_Matrix  => Real_Matrix,
 516                              Result_Matrix => Complex_Matrix,
 517                              Operation     => Compose_From_Cartesian);
 518
 519       ------------------------
 520       -- Compose_From_Polar --
 521       ------------------------
 522
 523       function Compose_From_Polar is
 524         new Vector_Vector_Elementwise_Operation
 525                             (Left_Scalar   => Real'Base,
 526                              Right_Scalar  => Real'Base,
 527                              Result_Scalar => Complex,
 528                              Left_Vector   => Real_Vector,
 529                              Right_Vector  => Real_Vector,
 530                              Result_Vector => Complex_Vector,
 531                              Operation     => Compose_From_Polar);
 532
 533       function Compose_From_Polar is
 534         new Vector_Vector_Scalar_Elementwise_Operation
 535                             (X_Scalar      => Real'Base,
 536                              Y_Scalar      => Real'Base,
 537                              Z_Scalar      => Real'Base,
 538                              Result_Scalar => Complex,
 539                              X_Vector      => Real_Vector,
 540                              Y_Vector      => Real_Vector,
 541                              Result_Vector => Complex_Vector,
 542                              Operation     => Compose_From_Polar);
 543
 544       function Compose_From_Polar is
 545         new Matrix_Matrix_Elementwise_Operation
 546                             (Left_Scalar   => Real'Base,
 547                              Right_Scalar  => Real'Base,
 548                              Result_Scalar => Complex,
 549                              Left_Matrix   => Real_Matrix,
 550                              Right_Matrix  => Real_Matrix,
 551                              Result_Matrix => Complex_Matrix,
 552                              Operation     => Compose_From_Polar);
 553
 554       function Compose_From_Polar is
 555         new Matrix_Matrix_Scalar_Elementwise_Operation
 556                             (X_Scalar      => Real'Base,
 557                              Y_Scalar      => Real'Base,
 558                              Z_Scalar      => Real'Base,
 559                              Result_Scalar => Complex,
 560                              X_Matrix      => Real_Matrix,
 561                              Y_Matrix      => Real_Matrix,
 562                              Result_Matrix => Complex_Matrix,
 563                              Operation     => Compose_From_Polar);
 564
 565       ---------------
 566       -- Conjugate --
 567       ---------------
 568
 569       function Conjugate is new Vector_Elementwise_Operation
 570                             (X_Scalar      => Complex,
 571                              Result_Scalar => Complex,
 572                              X_Vector      => Complex_Vector,
 573                              Result_Vector => Complex_Vector,
 574                              Operation     => Conjugate);
 575
 576       function Conjugate is new Matrix_Elementwise_Operation
 577                             (X_Scalar      => Complex,
 578                              Result_Scalar => Complex,
 579                              X_Matrix      => Complex_Matrix,
 580                              Result_Matrix => Complex_Matrix,
 581                              Operation     => Conjugate);
 582
 583       --------
 584       -- Im --
 585       --------
 586
 587       function Im is new Vector_Elementwise_Operation
 588                             (X_Scalar      => Complex,
 589                              Result_Scalar => Real'Base,
 590                              X_Vector      => Complex_Vector,
 591                              Result_Vector => Real_Vector,
 592                              Operation     => Im);
 593
 594       function Im is new Matrix_Elementwise_Operation
 595                             (X_Scalar      => Complex,
 596                              Result_Scalar => Real'Base,
 597                              X_Matrix      => Complex_Matrix,
 598                              Result_Matrix => Real_Matrix,
 599                              Operation     => Im);
 600
 601       -------------
 602       -- Modulus --
 603       -------------
 604
 605       function Modulus is new Vector_Elementwise_Operation
 606                             (X_Scalar      => Complex,
 607                              Result_Scalar => Real'Base,
 608                              X_Vector      => Complex_Vector,
 609                              Result_Vector => Real_Vector,
 610                              Operation     => Modulus);
 611
 612       function Modulus is new Matrix_Elementwise_Operation
 613                             (X_Scalar      => Complex,
 614                              Result_Scalar => Real'Base,
 615                              X_Matrix      => Complex_Matrix,
 616                              Result_Matrix => Real_Matrix,
 617                              Operation     => Modulus);
 618
 619       --------
 620       -- Re --
 621       --------
 622
 623       function Re is new Vector_Elementwise_Operation
 624                             (X_Scalar      => Complex,
 625                              Result_Scalar => Real'Base,
 626                              X_Vector      => Complex_Vector,
 627                              Result_Vector => Real_Vector,
 628                              Operation     => Re);
 629
 630       function Re is new Matrix_Elementwise_Operation
 631                             (X_Scalar      => Complex,
 632                              Result_Scalar => Real'Base,
 633                              X_Matrix      => Complex_Matrix,
 634                              Result_Matrix => Real_Matrix,
 635                              Operation     => Re);
 636
 637       ------------
 638       -- Set_Im --
 639       ------------
 640
 641       procedure Set_Im is new Update_Vector_With_Vector
 642                             (X_Scalar      => Complex,
 643                              Y_Scalar      => Real'Base,
 644                              X_Vector      => Complex_Vector,
 645                              Y_Vector      => Real_Vector,
 646                              Update        => Set_Im);
 647
 648       procedure Set_Im is new Update_Matrix_With_Matrix
 649                             (X_Scalar      => Complex,
 650                              Y_Scalar      => Real'Base,
 651                              X_Matrix      => Complex_Matrix,
 652                              Y_Matrix      => Real_Matrix,
 653                              Update        => Set_Im);
 654
 655       ------------
 656       -- Set_Re --
 657       ------------
 658
 659       procedure Set_Re is new Update_Vector_With_Vector
 660                             (X_Scalar      => Complex,
 661                              Y_Scalar      => Real'Base,
 662                              X_Vector      => Complex_Vector,
 663                              Y_Vector      => Real_Vector,
 664                              Update        => Set_Re);
 665
 666       procedure Set_Re is new Update_Matrix_With_Matrix
 667                             (X_Scalar      => Complex,
 668                              Y_Scalar      => Real'Base,
 669                              X_Matrix      => Complex_Matrix,
 670                              Y_Matrix      => Real_Matrix,
 671                              Update        => Set_Re);
 672
 673       -----------------
 674       -- Unit_Matrix --
 675       -----------------
 676
 677       function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix
 678                             (Scalar        => Complex,
 679                              Matrix        => Complex_Matrix,
 680                              Zero          => (0.0, 0.0),
 681                              One           => (1.0, 0.0));
 682
 683       function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector
 684                             (Scalar        => Complex,
 685                              Vector        => Complex_Vector,
 686                              Zero          => (0.0, 0.0),
 687                              One           => (1.0, 0.0));
 688
 689    end Instantiations;
 690
 691    ---------
 692    -- "*" --
 693    ---------
 694
 695    function "*"
 696      (Left  : Complex_Vector;
 697       Right : Complex_Vector) return Complex
 698    is
 699    begin
 700       if Left'Length /= Right'Length then
 701          raise Constraint_Error with
 702             "vectors are of different length in inner product";
 703       end if;
 704
 705       return dot (Left'Length, X => Left, Y => Right);
 706    end "*";
 707
 708    function "*"
 709      (Left  : Real_Vector;
 710       Right : Complex_Vector) return Complex
 711      renames Instantiations."*";
 712
 713    function "*"
 714      (Left  : Complex_Vector;
 715       Right : Real_Vector) return Complex
 716      renames Instantiations."*";
 717
 718    function "*"
 719      (Left  : Complex;
 720       Right : Complex_Vector) return Complex_Vector
 721      renames Instantiations."*";
 722
 723    function "*"
 724      (Left  : Complex_Vector;
 725       Right : Complex) return Complex_Vector
 726      renames Instantiations."*";
 727
 728    function "*"
 729      (Left  : Real'Base;
 730       Right : Complex_Vector) return Complex_Vector
 731      renames Instantiations."*";
 732
 733    function "*"
 734      (Left  : Complex_Vector;
 735       Right : Real'Base) return Complex_Vector
 736      renames Instantiations."*";
 737
 738    function "*"
 739      (Left  : Complex_Matrix;
 740       Right : Complex_Matrix)
 741       return  Complex_Matrix
 742    is
 743       R : Complex_Matrix (Left'Range (1), Right'Range (2));
 744
 745    begin
 746       if Left'Length (2) /= Right'Length (1) then
 747          raise Constraint_Error with
 748             "incompatible dimensions in matrix-matrix multipication";
 749       end if;
 750
 751       gemm (Trans_A => No_Trans'Access,
 752             Trans_B => No_Trans'Access,
 753             M       => Right'Length (2),
 754             N       => Left'Length (1),
 755             K       => Right'Length (1),
 756             A       => Right,
 757             Ld_A    => Right'Length (2),
 758             B       => Left,
 759             Ld_B    => Left'Length (2),
 760             C       => R,
 761             Ld_C    => R'Length (2));
 762
 763       return R;
 764    end "*";
 765
 766    function "*"
 767      (Left  : Complex_Vector;
 768       Right : Complex_Vector) return Complex_Matrix
 769      renames Instantiations."*";
 770
 771    function "*"
 772      (Left  : Complex_Vector;
 773       Right : Complex_Matrix) return Complex_Vector
 774    is
 775       R : Complex_Vector (Right'Range (2));
 776
 777    begin
 778       if Left'Length /= Right'Length (1) then
 779          raise Constraint_Error with
 780            "incompatible dimensions in vector-matrix multiplication";
 781       end if;
 782
 783       gemv (Trans => No_Trans'Access,
 784             M     => Right'Length (2),
 785             N     => Right'Length (1),
 786             A     => Right,
 787             Ld_A  => Right'Length (2),
 788             X     => Left,
 789             Y     => R);
 790
 791       return R;
 792    end "*";
 793
 794    function "*"
 795      (Left  : Complex_Matrix;
 796       Right : Complex_Vector) return Complex_Vector
 797    is
 798       R : Complex_Vector (Left'Range (1));
 799
 800    begin
 801       if Left'Length (2) /= Right'Length then
 802          raise Constraint_Error with
 803             "incompatible dimensions in matrix-vector multiplication";
 804       end if;
 805
 806       gemv (Trans => Trans'Access,
 807             M     => Left'Length (2),
 808             N     => Left'Length (1),
 809             A     => Left,
 810             Ld_A  => Left'Length (2),
 811             X     => Right,
 812             Y     => R);
 813
 814       return R;
 815    end "*";
 816
 817    function "*"
 818      (Left  : Real_Matrix;
 819       Right : Complex_Matrix) return Complex_Matrix
 820      renames Instantiations."*";
 821
 822    function "*"
 823      (Left  : Complex_Matrix;
 824       Right : Real_Matrix) return Complex_Matrix
 825      renames Instantiations."*";
 826
 827    function "*"
 828      (Left  : Real_Vector;
 829       Right : Complex_Vector) return Complex_Matrix
 830      renames Instantiations."*";
 831
 832    function "*"
 833      (Left  : Complex_Vector;
 834       Right : Real_Vector) return Complex_Matrix
 835      renames Instantiations."*";
 836
 837    function "*"
 838      (Left  : Real_Vector;
 839       Right : Complex_Matrix) return Complex_Vector
 840      renames Instantiations."*";
 841
 842    function "*"
 843      (Left  : Complex_Vector;
 844       Right : Real_Matrix) return Complex_Vector
 845      renames Instantiations."*";
 846
 847    function "*"
 848      (Left  : Real_Matrix;
 849       Right : Complex_Vector) return Complex_Vector
 850      renames Instantiations."*";
 851
 852    function "*"
 853      (Left  : Complex_Matrix;
 854       Right : Real_Vector) return Complex_Vector
 855      renames Instantiations."*";
 856
 857    function "*"
 858      (Left  : Complex;
 859       Right : Complex_Matrix) return Complex_Matrix
 860      renames Instantiations."*";
 861
 862    function "*"
 863      (Left  : Complex_Matrix;
 864       Right : Complex) return Complex_Matrix
 865      renames Instantiations."*";
 866
 867    function "*"
 868      (Left  : Real'Base;
 869       Right : Complex_Matrix) return Complex_Matrix
 870      renames Instantiations."*";
 871
 872    function "*"
 873      (Left  : Complex_Matrix;
 874       Right : Real'Base) return Complex_Matrix
 875      renames Instantiations."*";
 876
 877    ---------
 878    -- "+" --
 879    ---------
 880
 881    function "+" (Right : Complex_Vector) return Complex_Vector
 882      renames Instantiations."+";
 883
 884    function "+"
 885      (Left  : Complex_Vector;
 886       Right : Complex_Vector) return Complex_Vector
 887      renames Instantiations."+";
 888
 889    function "+"
 890      (Left  : Real_Vector;
 891       Right : Complex_Vector) return Complex_Vector
 892      renames Instantiations."+";
 893
 894    function "+"
 895      (Left  : Complex_Vector;
 896       Right : Real_Vector) return Complex_Vector
 897      renames Instantiations."+";
 898
 899    function "+" (Right : Complex_Matrix) return Complex_Matrix
 900      renames Instantiations."+";
 901
 902    function "+"
 903      (Left  : Complex_Matrix;
 904       Right : Complex_Matrix) return Complex_Matrix
 905      renames Instantiations."+";
 906
 907    function "+"
 908      (Left  : Real_Matrix;
 909       Right : Complex_Matrix) return Complex_Matrix
 910      renames Instantiations."+";
 911
 912    function "+"
 913      (Left  : Complex_Matrix;
 914       Right : Real_Matrix) return Complex_Matrix
 915      renames Instantiations."+";
 916
 917    ---------
 918    -- "-" --
 919    ---------
 920
 921    function "-"
 922      (Right : Complex_Vector) return Complex_Vector
 923      renames Instantiations."-";
 924
 925    function "-"
 926      (Left  : Complex_Vector;
 927       Right : Complex_Vector) return Complex_Vector
 928      renames Instantiations."-";
 929
 930    function "-"
 931      (Left  : Real_Vector;
 932       Right : Complex_Vector) return Complex_Vector
 933       renames Instantiations."-";
 934
 935    function "-"
 936      (Left  : Complex_Vector;
 937       Right : Real_Vector) return Complex_Vector
 938      renames Instantiations."-";
 939
 940    function "-" (Right : Complex_Matrix) return Complex_Matrix
 941      renames Instantiations."-";
 942
 943    function "-"
 944      (Left  : Complex_Matrix;
 945       Right : Complex_Matrix) return Complex_Matrix
 946      renames Instantiations."-";
 947
 948    function "-"
 949      (Left  : Real_Matrix;
 950       Right : Complex_Matrix) return Complex_Matrix
 951      renames Instantiations."-";
 952
 953    function "-"
 954      (Left  : Complex_Matrix;
 955       Right : Real_Matrix) return Complex_Matrix
 956      renames Instantiations."-";
 957
 958    ---------
 959    -- "/" --
 960    ---------
 961
 962    function "/"
 963      (Left  : Complex_Vector;
 964       Right : Complex) return Complex_Vector
 965      renames Instantiations."/";
 966
 967    function "/"
 968      (Left  : Complex_Vector;
 969       Right : Real'Base) return Complex_Vector
 970      renames Instantiations."/";
 971
 972    function "/"
 973      (Left  : Complex_Matrix;
 974       Right : Complex) return Complex_Matrix
 975      renames Instantiations."/";
 976
 977    function "/"
 978      (Left  : Complex_Matrix;
 979       Right : Real'Base) return Complex_Matrix
 980      renames Instantiations."/";
 981
 982    -----------
 983    -- "abs" --
 984    -----------
 985
 986    function "abs" (Right : Complex_Vector) return Complex is
 987    begin
 988       return (nrm2 (Right'Length, Right), 0.0);
 989    end "abs";
 990
 991    --------------
 992    -- Argument --
 993    --------------
 994
 995    function Argument (X : Complex_Vector) return Real_Vector
 996      renames Instantiations.Argument;
 997
 998    function Argument
 999      (X     : Complex_Vector;
1000       Cycle : Real'Base) return Real_Vector
1001      renames Instantiations.Argument;
1002
1003    function Argument (X : Complex_Matrix) return Real_Matrix
1004      renames Instantiations.Argument;
1005
1006    function Argument
1007      (X     : Complex_Matrix;
1008       Cycle : Real'Base) return Real_Matrix
1009      renames Instantiations.Argument;
1010
1011    ----------------------------
1012    -- Compose_From_Cartesian --
1013    ----------------------------
1014
1015    function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector
1016      renames Instantiations.Compose_From_Cartesian;
1017
1018    function Compose_From_Cartesian
1019      (Re : Real_Vector;
1020       Im : Real_Vector) return Complex_Vector
1021      renames Instantiations.Compose_From_Cartesian;
1022
1023    function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix
1024      renames Instantiations.Compose_From_Cartesian;
1025
1026    function Compose_From_Cartesian
1027      (Re : Real_Matrix;
1028       Im : Real_Matrix) return Complex_Matrix
1029      renames Instantiations.Compose_From_Cartesian;
1030
1031    ------------------------
1032    -- Compose_From_Polar --
1033    ------------------------
1034
1035    function Compose_From_Polar
1036      (Modulus  : Real_Vector;
1037       Argument : Real_Vector) return Complex_Vector
1038      renames Instantiations.Compose_From_Polar;
1039
1040    function Compose_From_Polar
1041      (Modulus  : Real_Vector;
1042       Argument : Real_Vector;
1043       Cycle    : Real'Base) return Complex_Vector
1044      renames Instantiations.Compose_From_Polar;
1045
1046    function Compose_From_Polar
1047      (Modulus  : Real_Matrix;
1048       Argument : Real_Matrix) return Complex_Matrix
1049      renames Instantiations.Compose_From_Polar;
1050
1051    function Compose_From_Polar
1052      (Modulus  : Real_Matrix;
1053       Argument : Real_Matrix;
1054       Cycle    : Real'Base) return Complex_Matrix
1055      renames Instantiations.Compose_From_Polar;
1056
1057    ---------------
1058    -- Conjugate --
1059    ---------------
1060
1061    function Conjugate (X : Complex_Vector) return Complex_Vector
1062      renames Instantiations.Conjugate;
1063
1064    function Conjugate (X : Complex_Matrix) return Complex_Matrix
1065      renames Instantiations.Conjugate;
1066
1067    -----------------
1068    -- Determinant --
1069    -----------------
1070
1071    function Determinant (A : Complex_Matrix) return Complex is
1072       N    : constant Integer := Length (A);
1073       LU   : Complex_Matrix (1 .. N, 1 .. N) := A;
1074       Piv  : Integer_Vector (1 .. N);
1075       Info : aliased Integer := -1;
1076       Neg  : Boolean;
1077       Det  : Complex;
1078
1079    begin
1080       if N = 0 then
1081          return (0.0, 0.0);
1082       end if;
1083
1084       getrf (N, N, LU, N, Piv, Info'Access);
1085
1086       if Info /= 0 then
1087          raise Constraint_Error with "ill-conditioned matrix";
1088       end if;
1089
1090       Det := LU (1, 1);
1091       Neg := Piv (1) /= 1;
1092
1093       for J in 2 .. N loop
1094          Det := Det * LU (J, J);
1095          Neg := Neg xor (Piv (J) /= J);
1096       end loop;
1097
1098       if Neg then
1099          return -Det;
1100
1101       else
1102          return Det;
1103       end if;
1104    end Determinant;
1105
1106    -----------------
1107    -- Eigensystem --
1108    -----------------
1109
1110    procedure Eigensystem
1111      (A       : in Complex_Matrix;
1112       Values  : out Real_Vector;
1113       Vectors : out Complex_Matrix)
1114    is
1115       Job_Z    : aliased Character := 'V';
1116       Rng      : aliased Character := 'A';
1117       Uplo     : aliased Character := 'U';
1118
1119       N        : constant Natural := Length (A);
1120       W        : BLAS_Real_Vector (Values'Range);
1121       M        : Integer;
1122       B        : Complex_Matrix (1 .. N, 1 .. N);
1123       L_Work   : Complex_Vector (1 .. 1);
1124       LR_Work  : BLAS_Real_Vector (1 .. 1);
1125       LI_Work  : Integer_Vector (1 .. 1);
1126       I_Supp_Z : Integer_Vector (1 .. 2 * N);
1127       Info     : aliased Integer;
1128
1129    begin
1130       if Values'Length /= N then
1131          raise Constraint_Error with "wrong length for output vector";
1132       end if;
1133
1134       if Vectors'First (1) /= A'First (1)
1135         or else Vectors'Last (1) /= A'Last (1)
1136         or else Vectors'First (2) /= A'First (2)
1137         or else Vectors'Last (2) /= A'Last (2)
1138       then
1139          raise Constraint_Error with "wrong dimensions for output matrix";
1140       end if;
1141
1142       if N = 0 then
1143          return;
1144       end if;
1145
1146       --  Check for hermitian matrix ???
1147       --  Copy only required triangle ???
1148
1149       B := A;
1150
1151       --  Find size of work area
1152
1153       heevr
1154         (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1155          M        => M,
1156          W        => W,
1157          Z        => Vectors,
1158          Ld_Z     => N,
1159          I_Supp_Z => I_Supp_Z,
1160          Work     => L_Work,
1161          L_Work   => -1,
1162          R_Work   => LR_Work,
1163          LR_Work  => -1,
1164          I_Work   => LI_Work,
1165          LI_Work  => -1,
1166          Info     => Info'Access);
1167
1168       if Info /= 0 then
1169          raise Constraint_Error;
1170       end if;
1171
1172       declare
1173          Work   : Complex_Vector (1 .. Integer (L_Work (1).Re));
1174          R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
1175          I_Work : Integer_Vector (1 .. LI_Work (1));
1176
1177       begin
1178          heevr
1179            (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1180             M        => M,
1181             W        => W,
1182             Z        => Vectors,
1183             Ld_Z     => N,
1184             I_Supp_Z => I_Supp_Z,
1185             Work     => Work,
1186             L_Work   => Work'Length,
1187             R_Work   => R_Work,
1188             LR_Work  => LR_Work'Length,
1189             I_Work   => I_Work,
1190             LI_Work  => LI_Work'Length,
1191             Info     => Info'Access);
1192
1193          if Info /= 0 then
1194             raise Constraint_Error with "inverting non-Hermetian matrix";
1195          end if;
1196
1197          for J in Values'Range loop
1198             Values (J) := W (J);
1199          end loop;
1200       end;
1201    end Eigensystem;
1202
1203    -----------------
1204    -- Eigenvalues --
1205    -----------------
1206
1207    procedure Eigenvalues
1208      (A      : Complex_Matrix;
1209       Values : out Real_Vector)
1210    is
1211       Job_Z    : aliased Character := 'N';
1212       Rng      : aliased Character := 'A';
1213       Uplo     : aliased Character := 'U';
1214       N        : constant Natural := Length (A);
1215       B        : Complex_Matrix (1 .. N, 1 .. N) := A;
1216       Z        : Complex_Matrix (1 .. 1, 1 .. 1);
1217       W        : BLAS_Real_Vector (Values'Range);
1218       L_Work   : Complex_Vector (1 .. 1);
1219       LR_Work  : BLAS_Real_Vector (1 .. 1);
1220       LI_Work  : Integer_Vector (1 .. 1);
1221       I_Supp_Z : Integer_Vector (1 .. 2 * N);
1222       M        : Integer;
1223       Info     : aliased Integer;
1224
1225    begin
1226       if Values'Length /= N then
1227          raise Constraint_Error with "wrong length for output vector";
1228       end if;
1229
1230       if N = 0 then
1231          return;
1232       end if;
1233
1234       --  Check for hermitian matrix ???
1235
1236       --  Find size of work area
1237
1238       heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1239              M        => M,
1240              W        => W,
1241              Z        => Z,
1242              Ld_Z     => 1,
1243              I_Supp_Z => I_Supp_Z,
1244              Work     => L_Work,  L_Work  => -1,
1245              R_Work   => LR_Work, LR_Work => -1,
1246              I_Work   => LI_Work, LI_Work => -1,
1247              Info     => Info'Access);
1248
1249       if Info /= 0 then
1250          raise Constraint_Error;
1251       end if;
1252
1253       declare
1254          Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1255          R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
1256          I_Work : Integer_Vector (1 .. LI_Work (1));
1257       begin
1258          heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
1259                 M        => M,
1260                 W        => W,
1261                 Z        => Z,
1262                 Ld_Z     => 1,
1263                 I_Supp_Z => I_Supp_Z,
1264                 Work     => Work,   L_Work  => Work'Length,
1265                 R_Work   => R_Work, LR_Work => R_Work'Length,
1266                 I_Work   => I_Work, LI_Work => I_Work'Length,
1267                 Info     => Info'Access);
1268
1269          if Info /= 0 then
1270             raise Constraint_Error with "inverting singular matrix";
1271          end if;
1272
1273          for J in Values'Range loop
1274             Values (J) := W (J);
1275          end loop;
1276       end;
1277    end Eigenvalues;
1278
1279    function Eigenvalues (A : Complex_Matrix) return Real_Vector is
1280       R : Real_Vector (A'Range (1));
1281    begin
1282       Eigenvalues (A, R);
1283       return R;
1284    end Eigenvalues;
1285
1286    --------
1287    -- Im --
1288    --------
1289
1290    function Im (X : Complex_Vector) return Real_Vector
1291      renames Instantiations.Im;
1292
1293    function Im (X : Complex_Matrix) return Real_Matrix
1294      renames Instantiations.Im;
1295
1296    -------------
1297    -- Inverse --
1298    -------------
1299
1300    procedure Inverse (A : Complex_Matrix; R : out Complex_Matrix) is
1301       N      : constant Integer := Length (A);
1302       Piv    : Integer_Vector (1 .. N);
1303       L_Work : Complex_Vector (1 .. 1);
1304       Info   : aliased Integer := -1;
1305
1306    begin
1307       --  All computations are done using column-major order, but this works
1308       --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
1309
1310       R := A;
1311
1312       --  Compute LU decomposition
1313
1314       getrf (M      => N,
1315              N      => N,
1316              A      => R,
1317              Ld_A   => N,
1318              I_Piv  => Piv,
1319              Info   => Info'Access);
1320
1321       if Info /= 0 then
1322          raise Constraint_Error with "inverting singular matrix";
1323       end if;
1324
1325       --  Determine size of work area
1326
1327       getri (N      => N,
1328              A      => R,
1329              Ld_A   => N,
1330              I_Piv  => Piv,
1331              Work   => L_Work,
1332              L_Work => -1,
1333              Info   => Info'Access);
1334
1335       if Info /= 0 then
1336          raise Constraint_Error;
1337       end if;
1338
1339       declare
1340          Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
1341
1342       begin
1343          --  Compute inverse from LU decomposition
1344
1345          getri (N      => N,
1346                 A      => R,
1347                 Ld_A   => N,
1348                 I_Piv  => Piv,
1349                 Work   => Work,
1350                 L_Work => Work'Length,
1351                 Info   => Info'Access);
1352
1353          if Info /= 0 then
1354             raise Constraint_Error with "inverting singular matrix";
1355          end if;
1356
1357          --  ??? Should iterate with gerfs, based on implementation advice
1358       end;
1359    end Inverse;
1360
1361    function Inverse (A : Complex_Matrix) return Complex_Matrix is
1362       R : Complex_Matrix (A'Range (2), A'Range (1));
1363    begin
1364       Inverse (A, R);
1365       return R;
1366    end Inverse;
1367
1368    -------------
1369    -- Modulus --
1370    -------------
1371
1372    function Modulus (X : Complex_Vector) return Real_Vector
1373      renames Instantiations.Modulus;
1374
1375    function Modulus (X : Complex_Matrix) return Real_Matrix
1376      renames Instantiations.Modulus;
1377
1378    --------
1379    -- Re --
1380    --------
1381
1382    function Re (X : Complex_Vector) return Real_Vector
1383      renames Instantiations.Re;
1384
1385    function Re (X : Complex_Matrix) return Real_Matrix
1386      renames Instantiations.Re;
1387
1388    ------------
1389    -- Set_Im --
1390    ------------
1391
1392    procedure Set_Im
1393      (X  : in out Complex_Matrix;
1394       Im : Real_Matrix)
1395      renames Instantiations.Set_Im;
1396
1397    procedure Set_Im
1398      (X  : in out Complex_Vector;
1399       Im : Real_Vector)
1400      renames Instantiations.Set_Im;
1401
1402    ------------
1403    -- Set_Re --
1404    ------------
1405
1406    procedure Set_Re
1407      (X  : in out Complex_Matrix;
1408       Re : Real_Matrix)
1409      renames Instantiations.Set_Re;
1410
1411    procedure Set_Re
1412      (X  : in out Complex_Vector;
1413       Re : Real_Vector)
1414      renames Instantiations.Set_Re;
1415
1416    -----------
1417    -- Solve --
1418    -----------
1419
1420    procedure Solve
1421      (A : Complex_Matrix;
1422       X : Complex_Vector;
1423       B : out Complex_Vector)
1424    is
1425    begin
1426       if Length (A) /= X'Length then
1427          raise Constraint_Error with
1428            "incompatible matrix and vector dimensions";
1429       end if;
1430
1431       --  ??? Should solve directly, is faster and more accurate
1432
1433       B := Inverse (A) * X;
1434    end Solve;
1435
1436    procedure Solve
1437      (A : Complex_Matrix;
1438       X : Complex_Matrix;
1439       B : out Complex_Matrix)
1440    is
1441    begin
1442       if Length (A) /= X'Length (1) then
1443          raise Constraint_Error with "incompatible matrix dimensions";
1444       end if;
1445
1446       --  ??? Should solve directly, is faster and more accurate
1447
1448       B := Inverse (A) * X;
1449    end Solve;
1450
1451    function Solve
1452      (A : Complex_Matrix;
1453       X : Complex_Vector) return Complex_Vector
1454    is
1455       B : Complex_Vector (A'Range (2));
1456    begin
1457       Solve (A, X, B);
1458       return B;
1459    end Solve;
1460
1461    function Solve (A, X : Complex_Matrix) return Complex_Matrix is
1462       B : Complex_Matrix (A'Range (2), X'Range (2));
1463    begin
1464       Solve (A, X, B);
1465       return B;
1466    end Solve;
1467
1468    ---------------
1469    -- Transpose --
1470    ---------------
1471
1472    function Transpose
1473      (X : Complex_Matrix) return Complex_Matrix
1474    is
1475       R : Complex_Matrix (X'Range (2), X'Range (1));
1476    begin
1477       Transpose (X, R);
1478       return R;
1479    end Transpose;
1480
1481    -----------------
1482    -- Unit_Matrix --
1483    -----------------
1484
1485    function Unit_Matrix
1486      (Order   : Positive;
1487       First_1 : Integer := 1;
1488       First_2 : Integer := 1) return Complex_Matrix
1489      renames Instantiations.Unit_Matrix;
1490
1491    -----------------
1492    -- Unit_Vector --
1493    -----------------
1494
1495    function Unit_Vector
1496      (Index : Integer;
1497       Order : Positive;
1498       First : Integer := 1) return Complex_Vector
1499      renames Instantiations.Unit_Vector;
1500
1501 end Ada.Numerics.Generic_Complex_Arrays;