gcc/ada/s-gecobl.adb

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