Subversion Repositories openrisc_2011-10-31

------------------------------------------------------------------------------

--                                                                          --

--                         GNAT RUN-TIME COMPONENTS                         --

--                                                                          --

--                     ADA.NUMERICS.GENERIC_REAL_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; use System;

with System.Generic_Real_BLAS;

with System.Generic_Real_LAPACK;

with System.Generic_Array_Operations; use System.Generic_Array_Operations;

package body Ada.Numerics.Generic_Real_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.

   package BLAS is

      new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);

   package LAPACK is

      new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);

   use BLAS, LAPACK;

   --  Procedure versions of functions returning unconstrained values.

   --  This allows for inlining the function wrapper.

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

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

   procedure Solve     (A : Real_Matrix; X : Real_Vector; B : out Real_Vector);

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

   procedure Transpose is new

     Generic_Array_Operations.Transpose

       (Scalar        => Real'Base,

        Matrix        => Real_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 (Real'Base, Real_Matrix);

   --  Instantiating the following subprograms directly would lead to

   --  name clashes, so use a local package.

   package Instantiations is

      function "+" is new

        Vector_Elementwise_Operation

          (X_Scalar      => Real'Base,

           Result_Scalar => Real'Base,

           X_Vector      => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "+");

      function "+" is new

        Matrix_Elementwise_Operation

          (X_Scalar      => Real'Base,

           Result_Scalar => Real'Base,

           X_Matrix      => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "+");

      function "+" is new

        Vector_Vector_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Vector   => Real_Vector,

           Right_Vector  => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "+");

      function "+" is new

        Matrix_Matrix_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Matrix   => Real_Matrix,

           Right_Matrix  => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "+");

      function "-" is new

        Vector_Elementwise_Operation

          (X_Scalar      => Real'Base,

           Result_Scalar => Real'Base,

           X_Vector      => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "-");

      function "-" is new

        Matrix_Elementwise_Operation

          (X_Scalar      => Real'Base,

           Result_Scalar => Real'Base,

           X_Matrix      => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "-");

      function "-" is new

        Vector_Vector_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Vector   => Real_Vector,

           Right_Vector  => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "-");

      function "-" is new

        Matrix_Matrix_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Matrix   => Real_Matrix,

           Right_Matrix  => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "-");

      function "*" is new

        Scalar_Vector_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Right_Vector  => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "*");

      function "*" is new

        Scalar_Matrix_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Right_Matrix  => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "*");

      function "*" is new

        Vector_Scalar_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Vector   => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "*");

      function "*" is new

        Matrix_Scalar_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Matrix   => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "*");

      function "*" is new

        Outer_Product

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Vector   => Real_Vector,

           Right_Vector  => Real_Vector,

           Matrix        => Real_Matrix);

      function "/" is new

        Vector_Scalar_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Vector   => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "/");

      function "/" is new

        Matrix_Scalar_Elementwise_Operation

          (Left_Scalar   => Real'Base,

           Right_Scalar  => Real'Base,

           Result_Scalar => Real'Base,

           Left_Matrix   => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "/");

      function "abs" is new

        Vector_Elementwise_Operation

          (X_Scalar      => Real'Base,

           Result_Scalar => Real'Base,

           X_Vector      => Real_Vector,

           Result_Vector => Real_Vector,

           Operation     => "abs");

      function "abs" is new

        Matrix_Elementwise_Operation

          (X_Scalar      => Real'Base,

           Result_Scalar => Real'Base,

           X_Matrix      => Real_Matrix,

           Result_Matrix => Real_Matrix,

           Operation     => "abs");

      function Unit_Matrix is new

        Generic_Array_Operations.Unit_Matrix

          (Scalar        => Real'Base,

           Matrix        => Real_Matrix,

           Zero          => 0.0,

           One           => 1.0);

      function Unit_Vector is new

        Generic_Array_Operations.Unit_Vector

          (Scalar        => Real'Base,

           Vector        => Real_Vector,

           Zero          => 0.0,

           One           => 1.0);

   end Instantiations;

   ---------

   -- "+" --

   ---------

   function "+" (Right : Real_Vector) return Real_Vector

      renames Instantiations."+";

   function "+" (Right : Real_Matrix) return Real_Matrix

      renames Instantiations."+";

   function "+" (Left, Right : Real_Vector) return Real_Vector

      renames Instantiations."+";

   function "+" (Left, Right : Real_Matrix) return Real_Matrix

      renames Instantiations."+";

   ---------

   -- "-" --

   ---------

   function "-" (Right : Real_Vector) return Real_Vector

      renames Instantiations."-";

   function "-" (Right : Real_Matrix) return Real_Matrix

      renames Instantiations."-";

   function "-" (Left, Right : Real_Vector) return Real_Vector

      renames Instantiations."-";

   function "-" (Left, Right : Real_Matrix) return Real_Matrix

      renames Instantiations."-";

   ---------

   -- "*" --

   ---------

   --  Scalar multiplication

   function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector

      renames Instantiations."*";

   function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector

      renames Instantiations."*";

   function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix

      renames Instantiations."*";

   function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix

      renames Instantiations."*";

   --  Vector multiplication

   function "*" (Left, Right : Real_Vector) return Real'Base 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, Right : Real_Vector) return Real_Matrix

      renames Instantiations."*";

   function "*"

     (Left : Real_Vector;

      Right : Real_Matrix) return Real_Vector

   is

      R : Real_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 : Real_Matrix;

      Right : Real_Vector) return Real_Vector

   is

      R : Real_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 "*";

   --  Matrix Multiplication

   function "*" (Left, Right : Real_Matrix) return Real_Matrix is

      R : Real_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 : Real_Vector; Right : Real'Base) return Real_Vector

      renames Instantiations."/";

   function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix

      renames Instantiations."/";

   -----------

   -- "abs" --

   -----------

   function "abs" (Right : Real_Vector) return Real'Base is

   begin

      return nrm2 (Right'Length, Right);

   end "abs";

   function "abs" (Right : Real_Vector) return Real_Vector

      renames Instantiations."abs";

   function "abs" (Right : Real_Matrix) return Real_Matrix

      renames Instantiations."abs";

   -----------------

   -- Determinant --

   -----------------

   function Determinant (A : Real_Matrix) return Real'Base is

      N    : constant Integer := Length (A);

      LU   : Real_Matrix (1 .. N, 1 .. N) := A;

      Piv  : Integer_Vector (1 .. N);

      Info : aliased Integer := -1;

      Det  : Real := 1.0;

   begin

      getrf (M     => N,

             N     => N,

             A     => LU,

             Ld_A  => N,

             I_Piv => Piv,

             Info  => Info'Access);

      if Info /= 0 then

         raise Constraint_Error with "ill-conditioned matrix";

      end if;

      for J in 1 .. N loop

         Det := (if Piv (J) /= J then -Det * LU (J, J) else Det * LU (J, J));

      end loop;

      return Det;

   end Determinant;

   -----------------

   -- Eigensystem --

   -----------------

   procedure Eigensystem

     (A       : Real_Matrix;

      Values  : out Real_Vector;

      Vectors : out Real_Matrix)

   is

      N      : constant Natural := Length (A);

      Tau    : Real_Vector (1 .. N);

      L_Work : Real_Vector (1 .. 1);

      Info   : aliased Integer;

      E : Real_Vector (1 .. N);

      pragma Warnings (Off, E);

   begin

      if Values'Length /= N then

         raise Constraint_Error with "wrong length for output vector";

      end if;

      if N = 0 then

         return;

      end if;

      --  Initialize working matrix and check for symmetric input matrix

      Transpose (A, Vectors);

      if A /= Vectors then

         raise Argument_Error with "matrix not symmetric";

      end if;

      --  Compute size of additional working space

      sytrd (Uplo   => Lower'Access,

             N      => N,

             A      => Vectors,

             Ld_A   => N,

             D      => Values,

             E      => E,

             Tau    => Tau,

             Work   => L_Work,

             L_Work => -1,

             Info   => Info'Access);

      declare

         Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));

         pragma Warnings (Off, Work);

         Comp_Z : aliased constant Character := 'V';

      begin

         --  Reduce matrix to tridiagonal form

         sytrd (Uplo   => Lower'Access,

                N      => N,

                A      => Vectors,

                Ld_A   => A'Length (1),

                D      => Values,

                E      => E,

                Tau    => Tau,

                Work   => Work,

                L_Work => Work'Length,

                Info   => Info'Access);

         if Info /= 0 then

            raise Program_Error;

         end if;

         --  Generate the real orthogonal matrix determined by sytrd

         orgtr (Uplo   => Lower'Access,

                N      => N,

                A      => Vectors,

                Ld_A   => N,

                Tau    => Tau,

                Work   => Work,

                L_Work => Work'Length,

                Info   => Info'Access);

         if Info /= 0 then

            raise Program_Error;

         end if;

         --  Compute all eigenvalues and eigenvectors using QR algorithm

         steqr (Comp_Z => Comp_Z'Access,

                N      => N,

                D      => Values,

                E      => E,

                Z      => Vectors,

                Ld_Z   => N,

                Work   => Work,

                Info   => Info'Access);

         if Info /= 0 then

            raise Constraint_Error with

               "eigensystem computation failed to converge";

         end if;

      end;

   end Eigensystem;

   -----------------

   -- Eigenvalues --

   -----------------

   procedure Eigenvalues

     (A      : Real_Matrix;

      Values : out Real_Vector)

   is

      N      : constant Natural := Length (A);

      L_Work : Real_Vector (1 .. 1);

      Info   : aliased Integer;

      B   : Real_Matrix (1 .. N, 1 .. N);

      Tau : Real_Vector (1 .. N);

      E   : Real_Vector (1 .. N);

      pragma Warnings (Off, B);

      pragma Warnings (Off, Tau);

      pragma Warnings (Off, E);

   begin

      if Values'Length /= N then

         raise Constraint_Error with "wrong length for output vector";

      end if;

      if N = 0 then

         return;

      end if;

      --  Initialize working matrix and check for symmetric input matrix

      Transpose (A, B);

      if A /= B then

         raise Argument_Error with "matrix not symmetric";

      end if;

      --  Find size of work area

      sytrd (Uplo   => Lower'Access,

             N      => N,

             A      => B,

             Ld_A   => N,

             D      => Values,

             E      => E,

             Tau    => Tau,

             Work   => L_Work,

             L_Work => -1,

             Info   => Info'Access);

      declare

         Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));

         pragma Warnings (Off, Work);

      begin

         --  Reduce matrix to tridiagonal form

         sytrd (Uplo   => Lower'Access,

                N      => N,

                A      => B,

                Ld_A   => A'Length (1),

                D      => Values,

                E      => E,

                Tau    => Tau,

                Work   => Work,

                L_Work => Work'Length,

                Info   => Info'Access);

         if Info /= 0 then

            raise Constraint_Error;

         end if;

         --  Compute all eigenvalues using QR algorithm

         sterf (N      => N,

                D      => Values,

                E      => E,

                Info   => Info'Access);

         if Info /= 0 then

            raise Constraint_Error with

               "eigenvalues computation failed to converge";

         end if;

      end;

   end Eigenvalues;

   function Eigenvalues (A : Real_Matrix) return Real_Vector is

      R : Real_Vector (A'Range (1));

   begin

      Eigenvalues (A, R);

      return R;

   end Eigenvalues;

   -------------

   -- Inverse --

   -------------

   procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is

      N      : constant Integer := Length (A);

      Piv    : Integer_Vector (1 .. N);

      L_Work : Real_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 : Real_Vector (1 .. Integer (L_Work (1)));

         pragma Warnings (Off, Work);

      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 : Real_Matrix) return Real_Matrix is

      R : Real_Matrix (A'Range (2), A'Range (1));

   begin

      Inverse (A, R);

      return R;

   end Inverse;

   -----------

   -- Solve --

   -----------

   procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_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 : Real_Matrix; X : Real_Matrix; B : out Real_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 : Real_Matrix; X : Real_Vector) return Real_Vector is

      B : Real_Vector (A'Range (2));

   begin

      Solve (A, X, B);

      return B;

   end Solve;

   function Solve (A, X : Real_Matrix) return Real_Matrix is

      B : Real_Matrix (A'Range (2), X'Range (2));

   begin

      Solve (A, X, B);

      return B;

   end Solve;

   ---------------

   -- Transpose --

   ---------------

   function Transpose (X : Real_Matrix) return Real_Matrix is

      R : Real_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 Real_Matrix

     renames Instantiations.Unit_Matrix;

   -----------------

   -- Unit_Vector --

   -----------------

   function Unit_Vector

     (Index : Integer;

      Order : Positive;

      First : Integer := 1) return Real_Vector

     renames Instantiations.Unit_Vector;

end Ada.Numerics.Generic_Real_Arrays;