gcc/ada/s-gerebl.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                         SYSTEM.GENERIC_REAL_BLAS                         --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --            Copyright (C) 2006, Free Software Foundation, Inc.            --
  10 --                                                                          --
  11 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12 -- terms of the  GNU General Public License as published  by the Free Soft- --
  13 -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  14 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  15 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  16 -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  17 -- for  more details.  You should have  received  a copy of the GNU General --
  18 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  19 -- to  the  Free Software Foundation,  51  Franklin  Street,  Fifth  Floor, --
  20 -- Boston, MA 02110-1301, USA.                                              --
  21 --                                                                          --
  22 -- As a special exception,  if other files  instantiate  generics from this --
  23 -- unit, or you link  this unit with other files  to produce an executable, --
  24 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  25 -- covered  by the  GNU  General  Public  License.  This exception does not --
  26 -- however invalidate  any other reasons why  the executable file  might be --
  27 -- covered by the  GNU Public License.                                      --
  28 --                                                                          --
  29 -- GNAT was originally developed  by the GNAT team at  New York University. --
  30 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  31 --                                                                          --
  32 ------------------------------------------------------------------------------
  33
  34 with Ada.Unchecked_Conversion;        use Ada;
  35 with Interfaces;                      use Interfaces;
  36 with Interfaces.Fortran;              use Interfaces.Fortran;
  37 with Interfaces.Fortran.BLAS;         use Interfaces.Fortran.BLAS;
  38 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
  39
  40 package body System.Generic_Real_BLAS is
  41
  42    Is_Single : constant Boolean :=
  43                  Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa
  44                   and then Fortran.Real (Real'First) = Fortran.Real'First
  45                   and then Fortran.Real (Real'Last) = Fortran.Real'Last;
  46
  47    Is_Double : constant Boolean :=
  48                  Real'Machine_Mantissa = Double_Precision'Machine_Mantissa
  49                   and then
  50                     Double_Precision (Real'First) = Double_Precision'First
  51                   and then
  52                     Double_Precision (Real'Last) = Double_Precision'Last;
  53
  54    --  Local subprograms
  55
  56    function To_Double_Precision (X : Real) return Double_Precision;
  57    pragma Inline_Always (To_Double_Precision);
  58
  59    function To_Real (X : Double_Precision) return Real;
  60    pragma Inline_Always (To_Real);
  61
  62    --  Instantiations
  63
  64    function To_Double_Precision is new
  65      Vector_Elementwise_Operation
  66        (X_Scalar      => Real,
  67         Result_Scalar => Double_Precision,
  68         X_Vector      => Real_Vector,
  69         Result_Vector => Double_Precision_Vector,
  70         Operation     => To_Double_Precision);
  71
  72    function To_Real is new
  73      Vector_Elementwise_Operation
  74        (X_Scalar      => Double_Precision,
  75         Result_Scalar => Real,
  76         X_Vector      => Double_Precision_Vector,
  77         Result_Vector => Real_Vector,
  78         Operation     => To_Real);
  79
  80    function To_Double_Precision is new
  81      Matrix_Elementwise_Operation
  82        (X_Scalar      => Real,
  83         Result_Scalar => Double_Precision,
  84         X_Matrix      => Real_Matrix,
  85         Result_Matrix => Double_Precision_Matrix,
  86         Operation     => To_Double_Precision);
  87
  88    function To_Real is new
  89      Matrix_Elementwise_Operation
  90        (X_Scalar      => Double_Precision,
  91         Result_Scalar => Real,
  92         X_Matrix      => Double_Precision_Matrix,
  93         Result_Matrix => Real_Matrix,
  94         Operation     => To_Real);
  95
  96    function To_Double_Precision (X : Real) return Double_Precision is
  97    begin
  98       return Double_Precision (X);
  99    end To_Double_Precision;
 100
 101    function To_Real (X : Double_Precision) return Real is
 102    begin
 103       return Real (X);
 104    end To_Real;
 105
 106    ---------
 107    -- dot --
 108    ---------
 109
 110    function dot
 111      (N     : Positive;
 112       X     : Real_Vector;
 113       Inc_X : Integer := 1;
 114       Y     : Real_Vector;
 115       Inc_Y : Integer := 1) return Real
 116    is
 117    begin
 118       if Is_Single then
 119          declare
 120             type X_Ptr is access all BLAS.Real_Vector (X'Range);
 121             type Y_Ptr is access all BLAS.Real_Vector (Y'Range);
 122             function Conv_X is new Unchecked_Conversion (Address, X_Ptr);
 123             function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
 124          begin
 125             return Real (sdot (N, Conv_X (X'Address).all, Inc_X,
 126                                   Conv_Y (Y'Address).all, Inc_Y));
 127          end;
 128
 129       elsif Is_Double then
 130          declare
 131             type X_Ptr is access all BLAS.Double_Precision_Vector (X'Range);
 132             type Y_Ptr is access all BLAS.Double_Precision_Vector (Y'Range);
 133             function Conv_X is new Unchecked_Conversion (Address, X_Ptr);
 134             function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
 135          begin
 136             return Real (ddot (N, Conv_X (X'Address).all, Inc_X,
 137                                   Conv_Y (Y'Address).all, Inc_Y));
 138          end;
 139
 140       else
 141          return Real (ddot (N, To_Double_Precision (X), Inc_X,
 142                                To_Double_Precision (Y), Inc_Y));
 143       end if;
 144    end dot;
 145
 146    ----------
 147    -- gemm --
 148    ----------
 149
 150    procedure gemm
 151      (Trans_A : access constant Character;
 152       Trans_B : access constant Character;
 153       M       : Positive;
 154       N       : Positive;
 155       K       : Positive;
 156       Alpha   : Real := 1.0;
 157       A       : Real_Matrix;
 158       Ld_A    : Integer;
 159       B       : Real_Matrix;
 160       Ld_B    : Integer;
 161       Beta    : Real := 0.0;
 162       C       : in out Real_Matrix;
 163       Ld_C    : Integer)
 164    is
 165    begin
 166       if Is_Single then
 167          declare
 168             subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2));
 169             subtype B_Type is BLAS.Real_Matrix (B'Range (1), B'Range (2));
 170             type C_Ptr is
 171               access all BLAS.Real_Matrix (C'Range (1), C'Range (2));
 172             function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
 173             function Conv_B is new Unchecked_Conversion (Real_Matrix, B_Type);
 174             function Conv_C is new Unchecked_Conversion (Address, C_Ptr);
 175          begin
 176             sgemm (Trans_A, Trans_B, M, N, K, Fortran.Real (Alpha),
 177                    Conv_A (A), Ld_A, Conv_B (B), Ld_B, Fortran.Real (Beta),
 178                    Conv_C (C'Address).all, Ld_C);
 179          end;
 180
 181       elsif Is_Double then
 182          declare
 183             subtype A_Type is
 184                Double_Precision_Matrix (A'Range (1), A'Range (2));
 185             subtype B_Type is
 186                Double_Precision_Matrix (B'Range (1), B'Range (2));
 187             type C_Ptr is
 188               access all Double_Precision_Matrix (C'Range (1), C'Range (2));
 189             function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
 190             function Conv_B is new Unchecked_Conversion (Real_Matrix, B_Type);
 191             function Conv_C is new Unchecked_Conversion (Address, C_Ptr);
 192          begin
 193             dgemm (Trans_A, Trans_B, M, N, K, Double_Precision (Alpha),
 194                    Conv_A (A), Ld_A, Conv_B (B), Ld_B, Double_Precision (Beta),
 195                    Conv_C (C'Address).all, Ld_C);
 196          end;
 197
 198       else
 199          declare
 200             DP_C : Double_Precision_Matrix (C'Range (1), C'Range (2));
 201          begin
 202             if Beta /= 0.0 then
 203                DP_C := To_Double_Precision (C);
 204             end if;
 205
 206             dgemm (Trans_A, Trans_B, M, N, K, Double_Precision (Alpha),
 207                    To_Double_Precision (A), Ld_A,
 208                    To_Double_Precision (B), Ld_B, Double_Precision (Beta),
 209                    DP_C, Ld_C);
 210
 211             C := To_Real (DP_C);
 212          end;
 213       end if;
 214    end gemm;
 215
 216    ----------
 217    -- gemv --
 218    ----------
 219
 220    procedure gemv
 221      (Trans : access constant Character;
 222       M     : Natural := 0;
 223       N     : Natural := 0;
 224       Alpha : Real := 1.0;
 225       A     : Real_Matrix;
 226       Ld_A  : Positive;
 227       X     : Real_Vector;
 228       Inc_X : Integer := 1;
 229       Beta  : Real := 0.0;
 230       Y     : in out Real_Vector;
 231       Inc_Y : Integer := 1)
 232    is
 233    begin
 234       if Is_Single then
 235          declare
 236             subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2));
 237             subtype X_Type is BLAS.Real_Vector (X'Range);
 238             type Y_Ptr is access all BLAS.Real_Vector (Y'Range);
 239             function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
 240             function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type);
 241             function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
 242          begin
 243             sgemv (Trans, M, N, Fortran.Real (Alpha),
 244                    Conv_A (A), Ld_A, Conv_X (X), Inc_X, Fortran.Real (Beta),
 245                    Conv_Y (Y'Address).all, Inc_Y);
 246          end;
 247
 248       elsif Is_Double then
 249          declare
 250             subtype A_Type is
 251                Double_Precision_Matrix (A'Range (1), A'Range (2));
 252             subtype X_Type is Double_Precision_Vector (X'Range);
 253             type Y_Ptr is access all Double_Precision_Vector (Y'Range);
 254             function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
 255             function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type);
 256             function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
 257          begin
 258             dgemv (Trans, M, N, Double_Precision (Alpha),
 259                    Conv_A (A), Ld_A, Conv_X (X), Inc_X,
 260                    Double_Precision (Beta),
 261                    Conv_Y (Y'Address).all, Inc_Y);
 262          end;
 263
 264       else
 265          declare
 266             DP_Y : Double_Precision_Vector (Y'Range);
 267          begin
 268             if Beta /= 0.0 then
 269                DP_Y := To_Double_Precision (Y);
 270             end if;
 271
 272             dgemv (Trans, M, N, Double_Precision (Alpha),
 273                    To_Double_Precision (A), Ld_A,
 274                    To_Double_Precision (X), Inc_X, Double_Precision (Beta),
 275                    DP_Y, Inc_Y);
 276
 277             Y := To_Real (DP_Y);
 278          end;
 279       end if;
 280    end gemv;
 281
 282    ----------
 283    -- nrm2 --
 284    ----------
 285
 286    function nrm2
 287      (N     : Natural;
 288       X     : Real_Vector;
 289       Inc_X : Integer := 1) return Real
 290    is
 291    begin
 292       if Is_Single then
 293          declare
 294             subtype X_Type is BLAS.Real_Vector (X'Range);
 295             function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type);
 296          begin
 297             return Real (snrm2 (N, Conv_X (X), Inc_X));
 298          end;
 299
 300       elsif Is_Double then
 301          declare
 302             subtype X_Type is Double_Precision_Vector (X'Range);
 303             function Conv_X is new Unchecked_Conversion (Real_Vector, X_Type);
 304          begin
 305             return Real (dnrm2 (N, Conv_X (X), Inc_X));
 306          end;
 307
 308       else
 309          return Real (dnrm2 (N, To_Double_Precision (X), Inc_X));
 310       end if;
 311    end nrm2;
 312
 313 end System.Generic_Real_BLAS;