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