gcc/ada/s-gerebl.adb

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