------------------------------------------------------------------------------
--                                                                          --
--                         GNAT COMPILER COMPONENTS                         --
--                                                                          --
--                   ADA.NUMERICS.GENERIC_COMPLEX_ARRAYS                    --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--            Copyright (C) 2006-2009, Free Software Foundation, Inc.       --
--                                                                          --
-- GNAT is free software;  you can  redistribute it  and/or modify it under --
-- terms of the  GNU General Public License as published  by the Free Soft- --
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
--                                                                          --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception,   --
-- version 3.1, as published by the Free Software Foundation.               --
--                                                                          --
-- You should have received a copy of the GNU General Public License and    --
-- a copy of the GCC Runtime Library Exception along with this program;     --
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
-- <http://www.gnu.org/licenses/>.                                          --
--                                                                          --
-- GNAT was originally developed  by the GNAT team at  New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
--                                                                          --
------------------------------------------------------------------------------

with System.Generic_Array_Operations; use System.Generic_Array_Operations;
with System.Generic_Complex_BLAS;
with System.Generic_Complex_LAPACK;

package body Ada.Numerics.Generic_Complex_Arrays is

   --  Operations involving inner products use BLAS library implementations.
   --  This allows larger matrices and vectors to be computed efficiently,
   --  taking into account memory hierarchy issues and vector instructions
   --  that vary widely between machines.

   --  Operations that are defined in terms of operations on the type Real,
   --  such as addition, subtraction and scaling, are computed in the canonical
   --  way looping over all elements.

   --  Operations for solving linear systems and computing determinant,
   --  eigenvalues, eigensystem and inverse, are implemented using the
   --  LAPACK library.

   type BLAS_Real_Vector is array (Integer range <>) of Real;

   package BLAS is new System.Generic_Complex_BLAS
     (Real           => Real,
      Complex_Types  => Complex_Types,
      Complex_Vector => Complex_Vector,
      Complex_Matrix => Complex_Matrix);

   package LAPACK is new System.Generic_Complex_LAPACK
     (Real           => Real,
      Real_Vector    => BLAS_Real_Vector,
      Complex_Types  => Complex_Types,
      Complex_Vector => Complex_Vector,
      Complex_Matrix => Complex_Matrix);

   subtype Real is Real_Arrays.Real;
   --  Work around visibility bug ???

   use BLAS, LAPACK;

   --  Procedure versions of functions returning unconstrained values.
   --  This allows for inlining the function wrapper.

   procedure Eigenvalues
     (A      : Complex_Matrix;
      Values : out Real_Vector);

   procedure Inverse
     (A      : Complex_Matrix;
      R      : out Complex_Matrix);

   procedure Solve
     (A      : Complex_Matrix;
      X      : Complex_Vector;
      B      : out Complex_Vector);

   procedure Solve
     (A      : Complex_Matrix;
      X      : Complex_Matrix;
      B      : out Complex_Matrix);

   procedure Transpose is new System.Generic_Array_Operations.Transpose
                                (Scalar => Complex,
                                 Matrix => Complex_Matrix);

   --  Helper function that raises a Constraint_Error is the argument is
   --  not a square matrix, and otherwise returns its length.

   function Length is new Square_Matrix_Length (Complex, Complex_Matrix);

   --  Instantiating the following subprograms directly would lead to
   --  name clashes, so use a local package.

   package Instantiations is

      ---------
      -- "*" --
      ---------

      function "*" is new Vector_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "*");

      function "*" is new Vector_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "*");

      function "*" is new Scalar_Vector_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "*");

      function "*" is new Scalar_Vector_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "*");

      function "*" is new Inner_Product
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Real_Vector,
                             Zero          => (0.0, 0.0));

      function "*" is new Inner_Product
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Right_Vector  => Complex_Vector,
                             Zero          => (0.0, 0.0));

      function "*" is new Outer_Product
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Complex_Vector,
                             Matrix        => Complex_Matrix);

      function "*" is new Outer_Product
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Right_Vector  => Complex_Vector,
                             Matrix        => Complex_Matrix);

      function "*" is new Outer_Product
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Real_Vector,
                             Matrix        => Complex_Matrix);

      function "*" is new Matrix_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "*");

      function "*" is new Matrix_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "*");

      function "*" is new Scalar_Matrix_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "*");

      function "*" is new Scalar_Matrix_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "*");

      function "*" is new Matrix_Vector_Product
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Matrix        => Real_Matrix,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Zero          => (0.0, 0.0));

      function "*" is new Matrix_Vector_Product
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Matrix        => Complex_Matrix,
                             Right_Vector  => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Zero          => (0.0, 0.0));

      function "*" is new Vector_Matrix_Product
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Matrix        => Complex_Matrix,
                             Result_Vector => Complex_Vector,
                             Zero          => (0.0, 0.0));

      function "*" is new Vector_Matrix_Product
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Matrix        => Real_Matrix,
                             Result_Vector => Complex_Vector,
                             Zero          => (0.0, 0.0));

      function "*" is new Matrix_Matrix_Product
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Real_Matrix,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Zero          => (0.0, 0.0));

      function "*" is new Matrix_Matrix_Product
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Right_Matrix  => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Zero          => (0.0, 0.0));

      ---------
      -- "+" --
      ---------

      function "+" is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Complex,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "+");

      function "+" is new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "+");

      function "+" is new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "+");

      function "+" is new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "+");

      function "+" is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Complex,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "+");

      function "+" is new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "+");

      function "+" is new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Real_Matrix,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "+");

      function "+" is new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Right_Matrix  => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "+");

      ---------
      -- "-" --
      ---------

      function "-" is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Complex,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "-");

      function "-" is new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "-");

      function "-" is new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Right_Vector  => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "-");

      function "-" is new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Right_Vector  => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "-");

      function "-" is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Complex,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "-");

      function "-" is new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "-");

      function "-" is new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Real_Matrix,
                             Right_Matrix  => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "-");

      function "-" is new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Right_Matrix  => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "-");

      ---------
      -- "/" --
      ---------

      function "/" is new Vector_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "/");

      function "/" is new Vector_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => "/");

      function "/" is new Matrix_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Complex,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "/");

      function "/" is new Matrix_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => "/");

      --------------
      -- Argument --
      --------------

      function Argument is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Real_Vector,
                             Operation     => Argument);

      function Argument is new Vector_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Real'Base,
                             Left_Vector   => Complex_Vector,
                             Result_Vector => Real_Vector,
                             Operation     => Argument);

      function Argument is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Real_Matrix,
                             Operation     => Argument);

      function Argument is new Matrix_Scalar_Elementwise_Operation
                            (Left_Scalar   => Complex,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Real'Base,
                             Left_Matrix   => Complex_Matrix,
                             Result_Matrix => Real_Matrix,
                             Operation     => Argument);

      ----------------------------
      -- Compose_From_Cartesian --
      ----------------------------

      function Compose_From_Cartesian is new Vector_Elementwise_Operation
                            (X_Scalar      => Real'Base,
                             Result_Scalar => Complex,
                             X_Vector      => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => Compose_From_Cartesian);

      function Compose_From_Cartesian is
         new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Right_Vector  => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => Compose_From_Cartesian);

      function Compose_From_Cartesian is new Matrix_Elementwise_Operation
                            (X_Scalar      => Real'Base,
                             Result_Scalar => Complex,
                             X_Matrix      => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => Compose_From_Cartesian);

      function Compose_From_Cartesian is
         new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Real_Matrix,
                             Right_Matrix  => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => Compose_From_Cartesian);

      ------------------------
      -- Compose_From_Polar --
      ------------------------

      function Compose_From_Polar is
        new Vector_Vector_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Vector   => Real_Vector,
                             Right_Vector  => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => Compose_From_Polar);

      function Compose_From_Polar is
        new Vector_Vector_Scalar_Elementwise_Operation
                            (X_Scalar      => Real'Base,
                             Y_Scalar      => Real'Base,
                             Z_Scalar      => Real'Base,
                             Result_Scalar => Complex,
                             X_Vector      => Real_Vector,
                             Y_Vector      => Real_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => Compose_From_Polar);

      function Compose_From_Polar is
        new Matrix_Matrix_Elementwise_Operation
                            (Left_Scalar   => Real'Base,
                             Right_Scalar  => Real'Base,
                             Result_Scalar => Complex,
                             Left_Matrix   => Real_Matrix,
                             Right_Matrix  => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => Compose_From_Polar);

      function Compose_From_Polar is
        new Matrix_Matrix_Scalar_Elementwise_Operation
                            (X_Scalar      => Real'Base,
                             Y_Scalar      => Real'Base,
                             Z_Scalar      => Real'Base,
                             Result_Scalar => Complex,
                             X_Matrix      => Real_Matrix,
                             Y_Matrix      => Real_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => Compose_From_Polar);

      ---------------
      -- Conjugate --
      ---------------

      function Conjugate is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Complex,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Complex_Vector,
                             Operation     => Conjugate);

      function Conjugate is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Complex,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Complex_Matrix,
                             Operation     => Conjugate);

      --------
      -- Im --
      --------

      function Im is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Real_Vector,
                             Operation     => Im);

      function Im is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Real_Matrix,
                             Operation     => Im);

      -------------
      -- Modulus --
      -------------

      function Modulus is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Real_Vector,
                             Operation     => Modulus);

      function Modulus is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Real_Matrix,
                             Operation     => Modulus);

      --------
      -- Re --
      --------

      function Re is new Vector_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Vector      => Complex_Vector,
                             Result_Vector => Real_Vector,
                             Operation     => Re);

      function Re is new Matrix_Elementwise_Operation
                            (X_Scalar      => Complex,
                             Result_Scalar => Real'Base,
                             X_Matrix      => Complex_Matrix,
                             Result_Matrix => Real_Matrix,
                             Operation     => Re);

      ------------
      -- Set_Im --
      ------------

      procedure Set_Im is new Update_Vector_With_Vector
                            (X_Scalar      => Complex,
                             Y_Scalar      => Real'Base,
                             X_Vector      => Complex_Vector,
                             Y_Vector      => Real_Vector,
                             Update        => Set_Im);

      procedure Set_Im is new Update_Matrix_With_Matrix
                            (X_Scalar      => Complex,
                             Y_Scalar      => Real'Base,
                             X_Matrix      => Complex_Matrix,
                             Y_Matrix      => Real_Matrix,
                             Update        => Set_Im);

      ------------
      -- Set_Re --
      ------------

      procedure Set_Re is new Update_Vector_With_Vector
                            (X_Scalar      => Complex,
                             Y_Scalar      => Real'Base,
                             X_Vector      => Complex_Vector,
                             Y_Vector      => Real_Vector,
                             Update        => Set_Re);

      procedure Set_Re is new Update_Matrix_With_Matrix
                            (X_Scalar      => Complex,
                             Y_Scalar      => Real'Base,
                             X_Matrix      => Complex_Matrix,
                             Y_Matrix      => Real_Matrix,
                             Update        => Set_Re);

      -----------------
      -- Unit_Matrix --
      -----------------

      function Unit_Matrix is new System.Generic_Array_Operations.Unit_Matrix
                            (Scalar        => Complex,
                             Matrix        => Complex_Matrix,
                             Zero          => (0.0, 0.0),
                             One           => (1.0, 0.0));

      function Unit_Vector is new System.Generic_Array_Operations.Unit_Vector
                            (Scalar        => Complex,
                             Vector        => Complex_Vector,
                             Zero          => (0.0, 0.0),
                             One           => (1.0, 0.0));

   end Instantiations;

   ---------
   -- "*" --
   ---------

   function "*"
     (Left  : Complex_Vector;
      Right : Complex_Vector) return Complex
   is
   begin
      if Left'Length /= Right'Length then
         raise Constraint_Error with
            "vectors are of different length in inner product";
      end if;

      return dot (Left'Length, X => Left, Y => Right);
   end "*";

   function "*"
     (Left  : Real_Vector;
      Right : Complex_Vector) return Complex
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Vector;
      Right : Real_Vector) return Complex
     renames Instantiations."*";

   function "*"
     (Left  : Complex;
      Right : Complex_Vector) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Vector;
      Right : Complex) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Real'Base;
      Right : Complex_Vector) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Vector;
      Right : Real'Base) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Matrix;
      Right : Complex_Matrix)
      return  Complex_Matrix
   is
      R : Complex_Matrix (Left'Range (1), Right'Range (2));

   begin
      if Left'Length (2) /= Right'Length (1) then
         raise Constraint_Error with
            "incompatible dimensions in matrix-matrix multiplication";
      end if;

      gemm (Trans_A => No_Trans'Access,
            Trans_B => No_Trans'Access,
            M       => Right'Length (2),
            N       => Left'Length (1),
            K       => Right'Length (1),
            A       => Right,
            Ld_A    => Right'Length (2),
            B       => Left,
            Ld_B    => Left'Length (2),
            C       => R,
            Ld_C    => R'Length (2));

      return R;
   end "*";

   function "*"
     (Left  : Complex_Vector;
      Right : Complex_Vector) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Vector;
      Right : Complex_Matrix) return Complex_Vector
   is
      R : Complex_Vector (Right'Range (2));

   begin
      if Left'Length /= Right'Length (1) then
         raise Constraint_Error with
           "incompatible dimensions in vector-matrix multiplication";
      end if;

      gemv (Trans => No_Trans'Access,
            M     => Right'Length (2),
            N     => Right'Length (1),
            A     => Right,
            Ld_A  => Right'Length (2),
            X     => Left,
            Y     => R);

      return R;
   end "*";

   function "*"
     (Left  : Complex_Matrix;
      Right : Complex_Vector) return Complex_Vector
   is
      R : Complex_Vector (Left'Range (1));

   begin
      if Left'Length (2) /= Right'Length then
         raise Constraint_Error with
            "incompatible dimensions in matrix-vector multiplication";
      end if;

      gemv (Trans => Trans'Access,
            M     => Left'Length (2),
            N     => Left'Length (1),
            A     => Left,
            Ld_A  => Left'Length (2),
            X     => Right,
            Y     => R);

      return R;
   end "*";

   function "*"
     (Left  : Real_Matrix;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Matrix;
      Right : Real_Matrix) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Real_Vector;
      Right : Complex_Vector) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Vector;
      Right : Real_Vector) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Real_Vector;
      Right : Complex_Matrix) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Vector;
      Right : Real_Matrix) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Real_Matrix;
      Right : Complex_Vector) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Matrix;
      Right : Real_Vector) return Complex_Vector
     renames Instantiations."*";

   function "*"
     (Left  : Complex;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Matrix;
      Right : Complex) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Real'Base;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."*";

   function "*"
     (Left  : Complex_Matrix;
      Right : Real'Base) return Complex_Matrix
     renames Instantiations."*";

   ---------
   -- "+" --
   ---------

   function "+" (Right : Complex_Vector) return Complex_Vector
     renames Instantiations."+";

   function "+"
     (Left  : Complex_Vector;
      Right : Complex_Vector) return Complex_Vector
     renames Instantiations."+";

   function "+"
     (Left  : Real_Vector;
      Right : Complex_Vector) return Complex_Vector
     renames Instantiations."+";

   function "+"
     (Left  : Complex_Vector;
      Right : Real_Vector) return Complex_Vector
     renames Instantiations."+";

   function "+" (Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."+";

   function "+"
     (Left  : Complex_Matrix;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."+";

   function "+"
     (Left  : Real_Matrix;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."+";

   function "+"
     (Left  : Complex_Matrix;
      Right : Real_Matrix) return Complex_Matrix
     renames Instantiations."+";

   ---------
   -- "-" --
   ---------

   function "-"
     (Right : Complex_Vector) return Complex_Vector
     renames Instantiations."-";

   function "-"
     (Left  : Complex_Vector;
      Right : Complex_Vector) return Complex_Vector
     renames Instantiations."-";

   function "-"
     (Left  : Real_Vector;
      Right : Complex_Vector) return Complex_Vector
      renames Instantiations."-";

   function "-"
     (Left  : Complex_Vector;
      Right : Real_Vector) return Complex_Vector
     renames Instantiations."-";

   function "-" (Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."-";

   function "-"
     (Left  : Complex_Matrix;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."-";

   function "-"
     (Left  : Real_Matrix;
      Right : Complex_Matrix) return Complex_Matrix
     renames Instantiations."-";

   function "-"
     (Left  : Complex_Matrix;
      Right : Real_Matrix) return Complex_Matrix
     renames Instantiations."-";

   ---------
   -- "/" --
   ---------

   function "/"
     (Left  : Complex_Vector;
      Right : Complex) return Complex_Vector
     renames Instantiations."/";

   function "/"
     (Left  : Complex_Vector;
      Right : Real'Base) return Complex_Vector
     renames Instantiations."/";

   function "/"
     (Left  : Complex_Matrix;
      Right : Complex) return Complex_Matrix
     renames Instantiations."/";

   function "/"
     (Left  : Complex_Matrix;
      Right : Real'Base) return Complex_Matrix
     renames Instantiations."/";

   -----------
   -- "abs" --
   -----------

   function "abs" (Right : Complex_Vector) return Complex is
   begin
      return (nrm2 (Right'Length, Right), 0.0);
   end "abs";

   --------------
   -- Argument --
   --------------

   function Argument (X : Complex_Vector) return Real_Vector
     renames Instantiations.Argument;

   function Argument
     (X     : Complex_Vector;
      Cycle : Real'Base) return Real_Vector
     renames Instantiations.Argument;

   function Argument (X : Complex_Matrix) return Real_Matrix
     renames Instantiations.Argument;

   function Argument
     (X     : Complex_Matrix;
      Cycle : Real'Base) return Real_Matrix
     renames Instantiations.Argument;

   ----------------------------
   -- Compose_From_Cartesian --
   ----------------------------

   function Compose_From_Cartesian (Re : Real_Vector) return Complex_Vector
     renames Instantiations.Compose_From_Cartesian;

   function Compose_From_Cartesian
     (Re : Real_Vector;
      Im : Real_Vector) return Complex_Vector
     renames Instantiations.Compose_From_Cartesian;

   function Compose_From_Cartesian (Re : Real_Matrix) return Complex_Matrix
     renames Instantiations.Compose_From_Cartesian;

   function Compose_From_Cartesian
     (Re : Real_Matrix;
      Im : Real_Matrix) return Complex_Matrix
     renames Instantiations.Compose_From_Cartesian;

   ------------------------
   -- Compose_From_Polar --
   ------------------------

   function Compose_From_Polar
     (Modulus  : Real_Vector;
      Argument : Real_Vector) return Complex_Vector
     renames Instantiations.Compose_From_Polar;

   function Compose_From_Polar
     (Modulus  : Real_Vector;
      Argument : Real_Vector;
      Cycle    : Real'Base) return Complex_Vector
     renames Instantiations.Compose_From_Polar;

   function Compose_From_Polar
     (Modulus  : Real_Matrix;
      Argument : Real_Matrix) return Complex_Matrix
     renames Instantiations.Compose_From_Polar;

   function Compose_From_Polar
     (Modulus  : Real_Matrix;
      Argument : Real_Matrix;
      Cycle    : Real'Base) return Complex_Matrix
     renames Instantiations.Compose_From_Polar;

   ---------------
   -- Conjugate --
   ---------------

   function Conjugate (X : Complex_Vector) return Complex_Vector
     renames Instantiations.Conjugate;

   function Conjugate (X : Complex_Matrix) return Complex_Matrix
     renames Instantiations.Conjugate;

   -----------------
   -- Determinant --
   -----------------

   function Determinant (A : Complex_Matrix) return Complex is
      N    : constant Integer := Length (A);
      LU   : Complex_Matrix (1 .. N, 1 .. N) := A;
      Piv  : Integer_Vector (1 .. N);
      Info : aliased Integer := -1;
      Neg  : Boolean;
      Det  : Complex;

   begin
      if N = 0 then
         return (0.0, 0.0);
      end if;

      getrf (N, N, LU, N, Piv, Info'Access);

      if Info /= 0 then
         raise Constraint_Error with "ill-conditioned matrix";
      end if;

      Det := LU (1, 1);
      Neg := Piv (1) /= 1;

      for J in 2 .. N loop
         Det := Det * LU (J, J);
         Neg := Neg xor (Piv (J) /= J);
      end loop;

      if Neg then
         return -Det;

      else
         return Det;
      end if;
   end Determinant;

   -----------------
   -- Eigensystem --
   -----------------

   procedure Eigensystem
     (A       : Complex_Matrix;
      Values  : out Real_Vector;
      Vectors : out Complex_Matrix)
   is
      Job_Z    : aliased Character := 'V';
      Rng      : aliased Character := 'A';
      Uplo     : aliased Character := 'U';

      N        : constant Natural := Length (A);
      W        : BLAS_Real_Vector (Values'Range);
      M        : Integer;
      B        : Complex_Matrix (1 .. N, 1 .. N);
      L_Work   : Complex_Vector (1 .. 1);
      LR_Work  : BLAS_Real_Vector (1 .. 1);
      LI_Work  : Integer_Vector (1 .. 1);
      I_Supp_Z : Integer_Vector (1 .. 2 * N);
      Info     : aliased Integer;

   begin
      if Values'Length /= N then
         raise Constraint_Error with "wrong length for output vector";
      end if;

      if Vectors'First (1) /= A'First (1)
        or else Vectors'Last (1) /= A'Last (1)
        or else Vectors'First (2) /= A'First (2)
        or else Vectors'Last (2) /= A'Last (2)
      then
         raise Constraint_Error with "wrong dimensions for output matrix";
      end if;

      if N = 0 then
         return;
      end if;

      --  Check for hermitian matrix ???
      --  Copy only required triangle ???

      B := A;

      --  Find size of work area

      heevr
        (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
         M        => M,
         W        => W,
         Z        => Vectors,
         Ld_Z     => N,
         I_Supp_Z => I_Supp_Z,
         Work     => L_Work,
         L_Work   => -1,
         R_Work   => LR_Work,
         LR_Work  => -1,
         I_Work   => LI_Work,
         LI_Work  => -1,
         Info     => Info'Access);

      if Info /= 0 then
         raise Constraint_Error;
      end if;

      declare
         Work   : Complex_Vector (1 .. Integer (L_Work (1).Re));
         R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
         I_Work : Integer_Vector (1 .. LI_Work (1));

      begin
         heevr
           (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
            M        => M,
            W        => W,
            Z        => Vectors,
            Ld_Z     => N,
            I_Supp_Z => I_Supp_Z,
            Work     => Work,
            L_Work   => Work'Length,
            R_Work   => R_Work,
            LR_Work  => LR_Work'Length,
            I_Work   => I_Work,
            LI_Work  => LI_Work'Length,
            Info     => Info'Access);

         if Info /= 0 then
            raise Constraint_Error with "inverting non-Hermitian matrix";
         end if;

         for J in Values'Range loop
            Values (J) := W (J);
         end loop;
      end;
   end Eigensystem;

   -----------------
   -- Eigenvalues --
   -----------------

   procedure Eigenvalues
     (A      : Complex_Matrix;
      Values : out Real_Vector)
   is
      Job_Z    : aliased Character := 'N';
      Rng      : aliased Character := 'A';
      Uplo     : aliased Character := 'U';
      N        : constant Natural := Length (A);
      B        : Complex_Matrix (1 .. N, 1 .. N) := A;
      Z        : Complex_Matrix (1 .. 1, 1 .. 1);
      W        : BLAS_Real_Vector (Values'Range);
      L_Work   : Complex_Vector (1 .. 1);
      LR_Work  : BLAS_Real_Vector (1 .. 1);
      LI_Work  : Integer_Vector (1 .. 1);
      I_Supp_Z : Integer_Vector (1 .. 2 * N);
      M        : Integer;
      Info     : aliased Integer;

   begin
      if Values'Length /= N then
         raise Constraint_Error with "wrong length for output vector";
      end if;

      if N = 0 then
         return;
      end if;

      --  Check for hermitian matrix ???

      --  Find size of work area

      heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
             M        => M,
             W        => W,
             Z        => Z,
             Ld_Z     => 1,
             I_Supp_Z => I_Supp_Z,
             Work     => L_Work,  L_Work  => -1,
             R_Work   => LR_Work, LR_Work => -1,
             I_Work   => LI_Work, LI_Work => -1,
             Info     => Info'Access);

      if Info /= 0 then
         raise Constraint_Error;
      end if;

      declare
         Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
         R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
         I_Work : Integer_Vector (1 .. LI_Work (1));
      begin
         heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
                M        => M,
                W        => W,
                Z        => Z,
                Ld_Z     => 1,
                I_Supp_Z => I_Supp_Z,
                Work     => Work,   L_Work  => Work'Length,
                R_Work   => R_Work, LR_Work => R_Work'Length,
                I_Work   => I_Work, LI_Work => I_Work'Length,
                Info     => Info'Access);

         if Info /= 0 then
            raise Constraint_Error with "inverting singular matrix";
         end if;

         for J in Values'Range loop
            Values (J) := W (J);
         end loop;
      end;
   end Eigenvalues;

   function Eigenvalues (A : Complex_Matrix) return Real_Vector is
      R : Real_Vector (A'Range (1));
   begin
      Eigenvalues (A, R);
      return R;
   end Eigenvalues;

   --------
   -- Im --
   --------

   function Im (X : Complex_Vector) return Real_Vector
     renames Instantiations.Im;

   function Im (X : Complex_Matrix) return Real_Matrix
     renames Instantiations.Im;

   -------------
   -- Inverse --
   -------------

   procedure Inverse (A : Complex_Matrix; R : out Complex_Matrix) is
      N      : constant Integer := Length (A);
      Piv    : Integer_Vector (1 .. N);
      L_Work : Complex_Vector (1 .. 1);
      Info   : aliased Integer := -1;

   begin
      --  All computations are done using column-major order, but this works
      --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).

      R := A;

      --  Compute LU decomposition

      getrf (M      => N,
             N      => N,
             A      => R,
             Ld_A   => N,
             I_Piv  => Piv,
             Info   => Info'Access);

      if Info /= 0 then
         raise Constraint_Error with "inverting singular matrix";
      end if;

      --  Determine size of work area

      getri (N      => N,
             A      => R,
             Ld_A   => N,
             I_Piv  => Piv,
             Work   => L_Work,
             L_Work => -1,
             Info   => Info'Access);

      if Info /= 0 then
         raise Constraint_Error;
      end if;

      declare
         Work : Complex_Vector (1 .. Integer (L_Work (1).Re));

      begin
         --  Compute inverse from LU decomposition

         getri (N      => N,
                A      => R,
                Ld_A   => N,
                I_Piv  => Piv,
                Work   => Work,
                L_Work => Work'Length,
                Info   => Info'Access);

         if Info /= 0 then
            raise Constraint_Error with "inverting singular matrix";
         end if;

         --  ??? Should iterate with gerfs, based on implementation advice
      end;
   end Inverse;

   function Inverse (A : Complex_Matrix) return Complex_Matrix is
      R : Complex_Matrix (A'Range (2), A'Range (1));
   begin
      Inverse (A, R);
      return R;
   end Inverse;

   -------------
   -- Modulus --
   -------------

   function Modulus (X : Complex_Vector) return Real_Vector
     renames Instantiations.Modulus;

   function Modulus (X : Complex_Matrix) return Real_Matrix
     renames Instantiations.Modulus;

   --------
   -- Re --
   --------

   function Re (X : Complex_Vector) return Real_Vector
     renames Instantiations.Re;

   function Re (X : Complex_Matrix) return Real_Matrix
     renames Instantiations.Re;

   ------------
   -- Set_Im --
   ------------

   procedure Set_Im
     (X  : in out Complex_Matrix;
      Im : Real_Matrix)
     renames Instantiations.Set_Im;

   procedure Set_Im
     (X  : in out Complex_Vector;
      Im : Real_Vector)
     renames Instantiations.Set_Im;

   ------------
   -- Set_Re --
   ------------

   procedure Set_Re
     (X  : in out Complex_Matrix;
      Re : Real_Matrix)
     renames Instantiations.Set_Re;

   procedure Set_Re
     (X  : in out Complex_Vector;
      Re : Real_Vector)
     renames Instantiations.Set_Re;

   -----------
   -- Solve --
   -----------

   procedure Solve
     (A : Complex_Matrix;
      X : Complex_Vector;
      B : out Complex_Vector)
   is
   begin
      if Length (A) /= X'Length then
         raise Constraint_Error with
           "incompatible matrix and vector dimensions";
      end if;

      --  ??? Should solve directly, is faster and more accurate

      B := Inverse (A) * X;
   end Solve;

   procedure Solve
     (A : Complex_Matrix;
      X : Complex_Matrix;
      B : out Complex_Matrix)
   is
   begin
      if Length (A) /= X'Length (1) then
         raise Constraint_Error with "incompatible matrix dimensions";
      end if;

      --  ??? Should solve directly, is faster and more accurate

      B := Inverse (A) * X;
   end Solve;

   function Solve
     (A : Complex_Matrix;
      X : Complex_Vector) return Complex_Vector
   is
      B : Complex_Vector (A'Range (2));
   begin
      Solve (A, X, B);
      return B;
   end Solve;

   function Solve (A, X : Complex_Matrix) return Complex_Matrix is
      B : Complex_Matrix (A'Range (2), X'Range (2));
   begin
      Solve (A, X, B);
      return B;
   end Solve;

   ---------------
   -- Transpose --
   ---------------

   function Transpose
     (X : Complex_Matrix) return Complex_Matrix
   is
      R : Complex_Matrix (X'Range (2), X'Range (1));
   begin
      Transpose (X, R);
      return R;
   end Transpose;

   -----------------
   -- Unit_Matrix --
   -----------------

   function Unit_Matrix
     (Order   : Positive;
      First_1 : Integer := 1;
      First_2 : Integer := 1) return Complex_Matrix
     renames Instantiations.Unit_Matrix;

   -----------------
   -- Unit_Vector --
   -----------------

   function Unit_Vector
     (Index : Integer;
      Order : Positive;
      First : Integer := 1) return Complex_Vector
     renames Instantiations.Unit_Vector;

end Ada.Numerics.Generic_Complex_Arrays;
