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-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_Complex_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    subtype Complex is Complex_Types.Complex;
  53
  54    --  Local subprograms
  55
  56    function To_Double_Precision (X : Real) return Double_Precision;
  57    pragma Inline (To_Double_Precision);
  58
  59    function To_Double_Complex (X : Complex) return Double_Complex;
  60    pragma Inline (To_Double_Complex);
  61
  62    function To_Complex (X : Double_Complex) return Complex;
  63    function To_Complex (X : Fortran.Complex) return Complex;
  64    pragma Inline (To_Complex);
  65
  66    function To_Fortran (X : Complex) return Fortran.Complex;
  67    pragma Inline (To_Fortran);
  68
  69    --  Instantiations
  70
  71    function To_Double_Complex is new
  72       Vector_Elementwise_Operation
  73        (X_Scalar      => Complex_Types.Complex,
  74         Result_Scalar => Fortran.Double_Complex,
  75         X_Vector      => Complex_Vector,
  76         Result_Vector => BLAS.Double_Complex_Vector,
  77         Operation     => To_Double_Complex);
  78
  79    function To_Complex is new
  80       Vector_Elementwise_Operation
  81        (X_Scalar      => Fortran.Double_Complex,
  82         Result_Scalar => Complex,
  83         X_Vector      => BLAS.Double_Complex_Vector,
  84         Result_Vector => Complex_Vector,
  85         Operation     => To_Complex);
  86
  87    function To_Double_Complex is new
  88       Matrix_Elementwise_Operation
  89        (X_Scalar      => Complex,
  90         Result_Scalar => Double_Complex,
  91         X_Matrix      => Complex_Matrix,
  92         Result_Matrix => BLAS.Double_Complex_Matrix,
  93         Operation     => To_Double_Complex);
  94
  95    function To_Complex is new
  96      Matrix_Elementwise_Operation
  97        (X_Scalar      => Double_Complex,
  98         Result_Scalar => Complex,
  99         X_Matrix      => BLAS.Double_Complex_Matrix,
 100         Result_Matrix => Complex_Matrix,
 101         Operation     => To_Complex);
 102
 103    function To_Double_Precision (X : Real) return Double_Precision is
 104    begin
 105       return Double_Precision (X);
 106    end To_Double_Precision;
 107
 108    function To_Double_Complex (X : Complex) return Double_Complex is
 109    begin
 110       return (To_Double_Precision (X.Re), To_Double_Precision (X.Im));
 111    end To_Double_Complex;
 112
 113    function To_Complex (X : Double_Complex) return Complex is
 114    begin
 115       return (Real (X.Re), Real (X.Im));
 116    end To_Complex;
 117
 118    function To_Complex (X : Fortran.Complex) return Complex is
 119    begin
 120       return (Real (X.Re), Real (X.Im));
 121    end To_Complex;
 122
 123    function To_Fortran (X : Complex) return Fortran.Complex is
 124    begin
 125       return (Fortran.Real (X.Re), Fortran.Real (X.Im));
 126    end To_Fortran;
 127
 128    ---------
 129    -- dot --
 130    ---------
 131
 132    function dot
 133      (N     : Positive;
 134       X     : Complex_Vector;
 135       Inc_X : Integer := 1;
 136       Y     : Complex_Vector;
 137       Inc_Y : Integer := 1) return Complex
 138    is
 139    begin
 140       if Is_Single then
 141          declare
 142             type X_Ptr is access all BLAS.Complex_Vector (X'Range);
 143             type Y_Ptr is access all BLAS.Complex_Vector (Y'Range);
 144             function Conv_X is new Unchecked_Conversion (Address, X_Ptr);
 145             function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
 146          begin
 147             return To_Complex (BLAS.cdotu (N, Conv_X (X'Address).all, Inc_X,
 148                                   Conv_Y (Y'Address).all, Inc_Y));
 149          end;
 150
 151       elsif Is_Double then
 152          declare
 153             type X_Ptr is access all BLAS.Double_Complex_Vector (X'Range);
 154             type Y_Ptr is access all BLAS.Double_Complex_Vector (Y'Range);
 155             function Conv_X is new Unchecked_Conversion (Address, X_Ptr);
 156             function Conv_Y is new Unchecked_Conversion (Address, Y_Ptr);
 157          begin
 158             return To_Complex (BLAS.zdotu (N, Conv_X (X'Address).all, Inc_X,
 159                                      Conv_Y (Y'Address).all, Inc_Y));
 160          end;
 161
 162       else
 163          return To_Complex (BLAS.zdotu (N, To_Double_Complex (X), Inc_X,
 164                                   To_Double_Complex (Y), Inc_Y));
 165       end if;
 166    end dot;
 167
 168    ----------
 169    -- gemm --
 170    ----------
 171
 172    procedure gemm
 173      (Trans_A : access constant Character;
 174       Trans_B : access constant Character;
 175       M       : Positive;
 176       N       : Positive;
 177       K       : Positive;
 178       Alpha   : Complex := (1.0, 0.0);
 179       A       : Complex_Matrix;
 180       Ld_A    : Integer;
 181       B       : Complex_Matrix;
 182       Ld_B    : Integer;
 183       Beta    : Complex := (0.0, 0.0);
 184       C       : in out Complex_Matrix;
 185       Ld_C    : Integer)
 186    is
 187    begin
 188       if Is_Single then
 189          declare
 190             subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2));
 191             subtype B_Type is BLAS.Complex_Matrix (B'Range (1), B'Range (2));
 192             type C_Ptr is
 193               access all BLAS.Complex_Matrix (C'Range (1), C'Range (2));
 194             function Conv_A is
 195                new Unchecked_Conversion (Complex_Matrix, A_Type);
 196             function Conv_B is
 197                new Unchecked_Conversion (Complex_Matrix, B_Type);
 198             function Conv_C is
 199                new Unchecked_Conversion (Address, C_Ptr);
 200          begin
 201             BLAS.cgemm (Trans_A, Trans_B, M, N, K, To_Fortran (Alpha),
 202                    Conv_A (A), Ld_A, Conv_B (B), Ld_B, To_Fortran (Beta),
 203                    Conv_C (C'Address).all, Ld_C);
 204          end;
 205
 206       elsif Is_Double then
 207          declare
 208             subtype A_Type is
 209                BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2));
 210             subtype B_Type is
 211                BLAS.Double_Complex_Matrix (B'Range (1), B'Range (2));
 212             type C_Ptr is access all
 213                BLAS.Double_Complex_Matrix (C'Range (1), C'Range (2));
 214             function Conv_A is
 215                new Unchecked_Conversion (Complex_Matrix, A_Type);
 216             function Conv_B is
 217                new Unchecked_Conversion (Complex_Matrix, B_Type);
 218             function Conv_C is new Unchecked_Conversion (Address, C_Ptr);
 219          begin
 220             BLAS.zgemm (Trans_A, Trans_B, M, N, K, To_Double_Complex (Alpha),
 221                    Conv_A (A), Ld_A, Conv_B (B), Ld_B,
 222                    To_Double_Complex (Beta),
 223                    Conv_C (C'Address).all, Ld_C);
 224          end;
 225
 226       else
 227          declare
 228             DP_C : BLAS.Double_Complex_Matrix (C'Range (1), C'Range (2));
 229          begin
 230             if Beta.Re /= 0.0 or else Beta.Im /= 0.0 then
 231                DP_C := To_Double_Complex (C);
 232             end if;
 233
 234             BLAS.zgemm (Trans_A, Trans_B, M, N, K, To_Double_Complex (Alpha),
 235                    To_Double_Complex (A), Ld_A,
 236                    To_Double_Complex (B), Ld_B, To_Double_Complex (Beta),
 237                    DP_C, Ld_C);
 238
 239             C := To_Complex (DP_C);
 240          end;
 241       end if;
 242    end gemm;
 243
 244    ----------
 245    -- gemv --
 246    ----------
 247
 248    procedure gemv
 249      (Trans : access constant Character;
 250       M     : Natural := 0;
 251       N     : Natural := 0;
 252       Alpha : Complex := (1.0, 0.0);
 253       A     : Complex_Matrix;
 254       Ld_A  : Positive;
 255       X     : Complex_Vector;
 256       Inc_X : Integer := 1;
 257       Beta  : Complex := (0.0, 0.0);
 258       Y     : in out Complex_Vector;
 259       Inc_Y : Integer := 1)
 260    is
 261    begin
 262       if Is_Single then
 263          declare
 264             subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2));
 265             subtype X_Type is BLAS.Complex_Vector (X'Range);
 266             type Y_Ptr is access all BLAS.Complex_Vector (Y'Range);
 267             function Conv_A is
 268                new Unchecked_Conversion (Complex_Matrix, A_Type);
 269             function Conv_X is
 270                new Unchecked_Conversion (Complex_Vector, X_Type);
 271             function Conv_Y is
 272                new Unchecked_Conversion (Address, Y_Ptr);
 273          begin
 274             BLAS.cgemv (Trans, M, N, To_Fortran (Alpha),
 275                    Conv_A (A), Ld_A, Conv_X (X), Inc_X, To_Fortran (Beta),
 276                    Conv_Y (Y'Address).all, Inc_Y);
 277          end;
 278
 279       elsif Is_Double then
 280          declare
 281             subtype A_Type is
 282                BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2));
 283             subtype X_Type is
 284                BLAS.Double_Complex_Vector (X'Range);
 285             type Y_Ptr is access all BLAS.Double_Complex_Vector (Y'Range);
 286             function Conv_A is
 287                new Unchecked_Conversion (Complex_Matrix, A_Type);
 288             function Conv_X is
 289                new Unchecked_Conversion (Complex_Vector, X_Type);
 290             function Conv_Y is
 291                new Unchecked_Conversion (Address, Y_Ptr);
 292          begin
 293             BLAS.zgemv (Trans, M, N, To_Double_Complex (Alpha),
 294                    Conv_A (A), Ld_A, Conv_X (X), Inc_X,
 295                    To_Double_Complex (Beta),
 296                    Conv_Y (Y'Address).all, Inc_Y);
 297          end;
 298
 299       else
 300          declare
 301             DP_Y : BLAS.Double_Complex_Vector (Y'Range);
 302          begin
 303             if Beta.Re /= 0.0 or else Beta.Im /= 0.0 then
 304                DP_Y := To_Double_Complex (Y);
 305             end if;
 306
 307             BLAS.zgemv (Trans, M, N, To_Double_Complex (Alpha),
 308                    To_Double_Complex (A), Ld_A,
 309                    To_Double_Complex (X), Inc_X, To_Double_Complex (Beta),
 310                    DP_Y, Inc_Y);
 311
 312             Y := To_Complex (DP_Y);
 313          end;
 314       end if;
 315    end gemv;
 316
 317    ----------
 318    -- nrm2 --
 319    ----------
 320
 321    function nrm2
 322      (N     : Natural;
 323       X     : Complex_Vector;
 324       Inc_X : Integer := 1) return Real
 325    is
 326    begin
 327       if Is_Single then
 328          declare
 329             subtype X_Type is BLAS.Complex_Vector (X'Range);
 330             function Conv_X is
 331                new Unchecked_Conversion (Complex_Vector, X_Type);
 332          begin
 333             return Real (BLAS.scnrm2 (N, Conv_X (X), Inc_X));
 334          end;
 335
 336       elsif Is_Double then
 337          declare
 338             subtype X_Type is BLAS.Double_Complex_Vector (X'Range);
 339             function Conv_X is
 340                new Unchecked_Conversion (Complex_Vector, X_Type);
 341          begin
 342             return Real (BLAS.dznrm2 (N, Conv_X (X), Inc_X));
 343          end;
 344
 345       else
 346          return Real (BLAS.dznrm2 (N, To_Double_Complex (X), Inc_X));
 347       end if;
 348    end nrm2;
 349
 350 end System.Generic_Complex_BLAS;