gcc/ada/s-gecola.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 _ L A P A C K        --
   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 Interfaces.Fortran.LAPACK;       use Interfaces.Fortran.LAPACK;
  37 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
  38
  39 package body System.Generic_Complex_LAPACK is
  40
  41    Is_Single : constant Boolean :=
  42                  Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa
  43                   and then Fortran.Real (Real'First) = Fortran.Real'First
  44                   and then Fortran.Real (Real'Last) = Fortran.Real'Last;
  45
  46    Is_Double : constant Boolean :=
  47                  Real'Machine_Mantissa = Double_Precision'Machine_Mantissa
  48                   and then
  49                     Double_Precision (Real'First) = Double_Precision'First
  50                   and then
  51                     Double_Precision (Real'Last) = Double_Precision'Last;
  52
  53    subtype Complex is Complex_Types.Complex;
  54
  55    --  Local subprograms
  56
  57    function To_Double_Precision (X : Real) return Double_Precision;
  58    pragma Inline (To_Double_Precision);
  59
  60    function To_Real (X : Double_Precision) return Real;
  61    pragma Inline (To_Real);
  62
  63    function To_Double_Complex (X : Complex) return Double_Complex;
  64    pragma Inline (To_Double_Complex);
  65
  66    function To_Complex (X : Double_Complex) return Complex;
  67    pragma Inline (To_Complex);
  68
  69    --  Instantiations
  70
  71    function To_Double_Precision is new
  72       Vector_Elementwise_Operation
  73        (X_Scalar      => Real,
  74         Result_Scalar => Double_Precision,
  75         X_Vector      => Real_Vector,
  76         Result_Vector => Double_Precision_Vector,
  77         Operation     => To_Double_Precision);
  78
  79    function To_Real is new
  80       Vector_Elementwise_Operation
  81        (X_Scalar      => Double_Precision,
  82         Result_Scalar => Real,
  83         X_Vector      => Double_Precision_Vector,
  84         Result_Vector => Real_Vector,
  85         Operation     => To_Real);
  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 => 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      => 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_Real (X : Double_Precision) return Real is
 109    begin
 110       return Real (X);
 111    end To_Real;
 112
 113    function To_Double_Complex (X : Complex) return Double_Complex is
 114    begin
 115       return (To_Double_Precision (X.Re), To_Double_Precision (X.Im));
 116    end To_Double_Complex;
 117
 118    function To_Complex (X : Double_Complex) return Complex is
 119    begin
 120       return (Real (X.Re), Real (X.Im));
 121    end To_Complex;
 122
 123    -----------
 124    -- getrf --
 125    -----------
 126
 127    procedure getrf
 128      (M     : Natural;
 129       N     : Natural;
 130       A     : in out Complex_Matrix;
 131       Ld_A  : Positive;
 132       I_Piv : out Integer_Vector;
 133       Info  : access Integer)
 134    is
 135    begin
 136       if Is_Single then
 137          declare
 138             type A_Ptr is
 139                access all BLAS.Complex_Matrix (A'Range (1), A'Range (2));
 140             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 141          begin
 142             cgetrf (M, N, Conv_A (A'Address).all, Ld_A,
 143                     LAPACK.Integer_Vector (I_Piv), Info);
 144          end;
 145
 146       elsif Is_Double then
 147          declare
 148             type A_Ptr is
 149                access all Double_Complex_Matrix (A'Range (1), A'Range (2));
 150             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 151          begin
 152             zgetrf (M, N, Conv_A (A'Address).all, Ld_A,
 153                     LAPACK.Integer_Vector (I_Piv), Info);
 154          end;
 155
 156       else
 157          declare
 158             DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
 159          begin
 160             DP_A := To_Double_Complex (A);
 161             zgetrf (M, N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), Info);
 162             A := To_Complex (DP_A);
 163          end;
 164       end if;
 165    end getrf;
 166
 167    -----------
 168    -- getri --
 169    -----------
 170
 171    procedure getri
 172      (N      : Natural;
 173       A      : in out Complex_Matrix;
 174       Ld_A   : Positive;
 175       I_Piv  : Integer_Vector;
 176       Work   : in out Complex_Vector;
 177       L_Work : Integer;
 178       Info   : access Integer)
 179    is
 180    begin
 181       if Is_Single then
 182          declare
 183             type A_Ptr is
 184                access all BLAS.Complex_Matrix (A'Range (1), A'Range (2));
 185             type Work_Ptr is
 186                access all BLAS.Complex_Vector (Work'Range);
 187             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 188             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 189          begin
 190             cgetri (N, Conv_A (A'Address).all, Ld_A,
 191                     LAPACK.Integer_Vector (I_Piv),
 192                     Conv_Work (Work'Address).all, L_Work,
 193                     Info);
 194          end;
 195
 196       elsif Is_Double then
 197          declare
 198             type A_Ptr is
 199                access all Double_Complex_Matrix (A'Range (1), A'Range (2));
 200             type Work_Ptr is
 201                access all Double_Complex_Vector (Work'Range);
 202             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 203             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 204          begin
 205             zgetri (N, Conv_A (A'Address).all, Ld_A,
 206                     LAPACK.Integer_Vector (I_Piv),
 207                     Conv_Work (Work'Address).all, L_Work,
 208                     Info);
 209          end;
 210
 211       else
 212          declare
 213             DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
 214             DP_Work : Double_Complex_Vector (Work'Range);
 215          begin
 216             DP_A := To_Double_Complex (A);
 217             zgetri (N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv),
 218                     DP_Work, L_Work, Info);
 219             A := To_Complex (DP_A);
 220             Work (1) := To_Complex (DP_Work (1));
 221          end;
 222       end if;
 223    end getri;
 224
 225    -----------
 226    -- getrs --
 227    -----------
 228
 229    procedure getrs
 230      (Trans  : access constant Character;
 231       N      : Natural;
 232       N_Rhs  : Natural;
 233       A      : Complex_Matrix;
 234       Ld_A   : Positive;
 235       I_Piv  : Integer_Vector;
 236       B      : in out Complex_Matrix;
 237       Ld_B   : Positive;
 238       Info   : access Integer)
 239    is
 240    begin
 241       if Is_Single then
 242          declare
 243             subtype A_Type is BLAS.Complex_Matrix (A'Range (1), A'Range (2));
 244             type B_Ptr is
 245                access all BLAS.Complex_Matrix (B'Range (1), B'Range (2));
 246             function Conv_A is
 247                new Unchecked_Conversion (Complex_Matrix, A_Type);
 248             function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
 249          begin
 250             cgetrs (Trans, N, N_Rhs,
 251                     Conv_A (A), Ld_A,
 252                     LAPACK.Integer_Vector (I_Piv),
 253                     Conv_B (B'Address).all, Ld_B,
 254                     Info);
 255          end;
 256
 257       elsif Is_Double then
 258          declare
 259             subtype A_Type is
 260                Double_Complex_Matrix (A'Range (1), A'Range (2));
 261             type B_Ptr is
 262                access all Double_Complex_Matrix (B'Range (1), B'Range (2));
 263             function Conv_A is
 264                new Unchecked_Conversion (Complex_Matrix, A_Type);
 265             function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
 266          begin
 267             zgetrs (Trans, N, N_Rhs,
 268                     Conv_A (A), Ld_A,
 269                     LAPACK.Integer_Vector (I_Piv),
 270                     Conv_B (B'Address).all, Ld_B,
 271                     Info);
 272          end;
 273
 274       else
 275          declare
 276             DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
 277             DP_B : Double_Complex_Matrix (B'Range (1), B'Range (2));
 278          begin
 279             DP_A := To_Double_Complex (A);
 280             DP_B := To_Double_Complex (B);
 281             zgetrs (Trans, N, N_Rhs,
 282                     DP_A, Ld_A,
 283                     LAPACK.Integer_Vector (I_Piv),
 284                     DP_B, Ld_B,
 285                     Info);
 286             B := To_Complex (DP_B);
 287          end;
 288       end if;
 289    end getrs;
 290
 291    procedure heevr
 292      (Job_Z    : access constant Character;
 293       Rng      : access constant Character;
 294       Uplo     : access constant Character;
 295       N        : Natural;
 296       A        : in out Complex_Matrix;
 297       Ld_A     : Positive;
 298       Vl, Vu   : Real := 0.0;
 299       Il, Iu   : Integer := 1;
 300       Abs_Tol  : Real := 0.0;
 301       M        : out Integer;
 302       W        : out Real_Vector;
 303       Z        : out Complex_Matrix;
 304       Ld_Z     : Positive;
 305       I_Supp_Z : out Integer_Vector;
 306       Work     : out Complex_Vector;
 307       L_Work   : Integer;
 308       R_Work   : out Real_Vector;
 309       LR_Work  : Integer;
 310       I_Work   : out Integer_Vector;
 311       LI_Work  : Integer;
 312       Info     : access Integer)
 313    is
 314    begin
 315       if Is_Single then
 316          declare
 317             type A_Ptr is
 318                access all BLAS.Complex_Matrix (A'Range (1), A'Range (2));
 319             type W_Ptr is
 320                access all BLAS.Real_Vector (W'Range);
 321             type Z_Ptr is
 322                access all BLAS.Complex_Matrix (Z'Range (1), Z'Range (2));
 323             type Work_Ptr is access all  BLAS.Complex_Vector (Work'Range);
 324             type R_Work_Ptr is access all BLAS.Real_Vector (R_Work'Range);
 325
 326             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 327             function Conv_W is new Unchecked_Conversion (Address, W_Ptr);
 328             function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
 329             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 330             function Conv_R_Work is
 331                new Unchecked_Conversion (Address, R_Work_Ptr);
 332          begin
 333             cheevr (Job_Z, Rng, Uplo, N,
 334                     Conv_A (A'Address).all, Ld_A,
 335                     Fortran.Real (Vl), Fortran.Real (Vu),
 336                     Il, Iu, Fortran.Real (Abs_Tol), M,
 337                     Conv_W (W'Address).all,
 338                     Conv_Z (Z'Address).all, Ld_Z,
 339                     LAPACK.Integer_Vector (I_Supp_Z),
 340                     Conv_Work (Work'Address).all, L_Work,
 341                     Conv_R_Work (R_Work'Address).all, LR_Work,
 342                     LAPACK.Integer_Vector (I_Work), LI_Work, Info);
 343          end;
 344
 345       elsif Is_Double then
 346          declare
 347             type A_Ptr is
 348               access all BLAS.Double_Complex_Matrix (A'Range (1), A'Range (2));
 349             type W_Ptr is
 350               access all BLAS.Double_Precision_Vector (W'Range);
 351             type Z_Ptr is
 352               access all BLAS.Double_Complex_Matrix (Z'Range (1), Z'Range (2));
 353             type Work_Ptr is
 354                access all BLAS.Double_Complex_Vector (Work'Range);
 355             type R_Work_Ptr is
 356                access all BLAS.Double_Precision_Vector (R_Work'Range);
 357
 358             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 359             function Conv_W is new Unchecked_Conversion (Address, W_Ptr);
 360             function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
 361             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 362             function Conv_R_Work is
 363                new Unchecked_Conversion (Address, R_Work_Ptr);
 364          begin
 365             zheevr (Job_Z, Rng, Uplo, N,
 366                     Conv_A (A'Address).all, Ld_A,
 367                     Double_Precision (Vl), Double_Precision (Vu),
 368                     Il, Iu, Double_Precision (Abs_Tol), M,
 369                     Conv_W (W'Address).all,
 370                     Conv_Z (Z'Address).all, Ld_Z,
 371                     LAPACK.Integer_Vector (I_Supp_Z),
 372                     Conv_Work (Work'Address).all, L_Work,
 373                     Conv_R_Work (R_Work'Address).all, LR_Work,
 374                     LAPACK.Integer_Vector (I_Work), LI_Work, Info);
 375          end;
 376
 377       else
 378          declare
 379             DP_A : Double_Complex_Matrix (A'Range (1), A'Range (2));
 380             DP_W : Double_Precision_Vector (W'Range);
 381             DP_Z : Double_Complex_Matrix (Z'Range (1), Z'Range (2));
 382             DP_Work : Double_Complex_Vector (Work'Range);
 383             DP_R_Work : Double_Precision_Vector (R_Work'Range);
 384
 385          begin
 386             DP_A := To_Double_Complex (A);
 387
 388             zheevr (Job_Z, Rng, Uplo, N,
 389                     DP_A, Ld_A,
 390                     Double_Precision (Vl), Double_Precision (Vu),
 391                     Il, Iu, Double_Precision (Abs_Tol), M,
 392                     DP_W, DP_Z, Ld_Z,
 393                     LAPACK.Integer_Vector (I_Supp_Z),
 394                     DP_Work, L_Work,
 395                     DP_R_Work, LR_Work,
 396                     LAPACK.Integer_Vector (I_Work), LI_Work, Info);
 397
 398             A := To_Complex (DP_A);
 399             W := To_Real (DP_W);
 400             Z := To_Complex (DP_Z);
 401
 402             Work (1) := To_Complex (DP_Work (1));
 403             R_Work (1) := To_Real (DP_R_Work (1));
 404          end;
 405       end if;
 406    end heevr;
 407
 408    -----------
 409    -- steqr --
 410    -----------
 411
 412    procedure steqr
 413      (Comp_Z : access constant Character;
 414       N      : Natural;
 415       D      : in out Real_Vector;
 416       E      : in out Real_Vector;
 417       Z      : in out Complex_Matrix;
 418       Ld_Z   : Positive;
 419       Work   : out Real_Vector;
 420       Info   : access Integer)
 421    is
 422    begin
 423       if Is_Single then
 424          declare
 425             type D_Ptr is access all BLAS.Real_Vector (D'Range);
 426             type E_Ptr is access all BLAS.Real_Vector (E'Range);
 427             type Z_Ptr is
 428                access all BLAS.Complex_Matrix (Z'Range (1), Z'Range (2));
 429             type Work_Ptr is
 430                access all BLAS.Real_Vector (Work'Range);
 431             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 432             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 433             function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
 434             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 435          begin
 436             csteqr (Comp_Z, N,
 437                     Conv_D (D'Address).all,
 438                     Conv_E (E'Address).all,
 439                     Conv_Z (Z'Address).all,
 440                     Ld_Z,
 441                     Conv_Work (Work'Address).all,
 442                     Info);
 443          end;
 444
 445       elsif Is_Double then
 446          declare
 447             type D_Ptr is access all Double_Precision_Vector (D'Range);
 448             type E_Ptr is access all Double_Precision_Vector (E'Range);
 449             type Z_Ptr is
 450                access all Double_Complex_Matrix (Z'Range (1), Z'Range (2));
 451             type Work_Ptr is
 452                access all Double_Precision_Vector (Work'Range);
 453             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 454             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 455             function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
 456             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 457          begin
 458             zsteqr (Comp_Z, N,
 459                     Conv_D (D'Address).all,
 460                     Conv_E (E'Address).all,
 461                     Conv_Z (Z'Address).all,
 462                     Ld_Z,
 463                     Conv_Work (Work'Address).all,
 464                     Info);
 465          end;
 466
 467       else
 468          declare
 469             DP_D    : Double_Precision_Vector (D'Range);
 470             DP_E    : Double_Precision_Vector (E'Range);
 471             DP_Z    : Double_Complex_Matrix (Z'Range (1), Z'Range (2));
 472             DP_Work : Double_Precision_Vector (Work'Range);
 473          begin
 474             DP_D := To_Double_Precision (D);
 475             DP_E := To_Double_Precision (E);
 476
 477             if Comp_Z.all = 'V' then
 478                DP_Z := To_Double_Complex (Z);
 479             end if;
 480
 481             zsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info);
 482
 483             D := To_Real (DP_D);
 484             E := To_Real (DP_E);
 485
 486             if Comp_Z.all /= 'N' then
 487                Z := To_Complex (DP_Z);
 488             end if;
 489          end;
 490       end if;
 491    end steqr;
 492
 493 end System.Generic_Complex_LAPACK;