gcc/ada/s-gecola.adb

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