gcc/ada/s-gerela.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                        SYSTEM.GENERIC_REAL_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_Real_LAPACK is
  42
  43    Is_Real : 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_Precision : constant Boolean :=
  49                            Real'Machine_Mantissa =
  50                                             Double_Precision'Machine_Mantissa
  51                             and then
  52                               Double_Precision (Real'First) =
  53                                             Double_Precision'First
  54                             and then
  55                               Double_Precision (Real'Last) =
  56                                             Double_Precision'Last;
  57
  58    --  Local subprograms
  59
  60    function To_Double_Precision (X : Real) return Double_Precision;
  61    pragma Inline_Always (To_Double_Precision);
  62
  63    function To_Real (X : Double_Precision) return Real;
  64    pragma Inline_Always (To_Real);
  65
  66    --  Instantiations
  67
  68    function To_Double_Precision is new
  69      Vector_Elementwise_Operation
  70        (X_Scalar      => Real,
  71         Result_Scalar => Double_Precision,
  72         X_Vector      => Real_Vector,
  73         Result_Vector => Double_Precision_Vector,
  74         Operation     => To_Double_Precision);
  75
  76    function To_Real is new
  77      Vector_Elementwise_Operation
  78        (X_Scalar      => Double_Precision,
  79         Result_Scalar => Real,
  80         X_Vector      => Double_Precision_Vector,
  81         Result_Vector => Real_Vector,
  82         Operation     => To_Real);
  83
  84    function To_Double_Precision is new
  85      Matrix_Elementwise_Operation
  86        (X_Scalar      => Real,
  87         Result_Scalar => Double_Precision,
  88         X_Matrix      => Real_Matrix,
  89         Result_Matrix => Double_Precision_Matrix,
  90         Operation     => To_Double_Precision);
  91
  92    function To_Real is new
  93      Matrix_Elementwise_Operation
  94        (X_Scalar      => Double_Precision,
  95         Result_Scalar => Real,
  96         X_Matrix      => Double_Precision_Matrix,
  97         Result_Matrix => Real_Matrix,
  98         Operation     => To_Real);
  99
 100    function To_Double_Precision (X : Real) return Double_Precision is
 101    begin
 102       return Double_Precision (X);
 103    end To_Double_Precision;
 104
 105    function To_Real (X : Double_Precision) return Real is
 106    begin
 107       return Real (X);
 108    end To_Real;
 109
 110    -----------
 111    -- getrf --
 112    -----------
 113
 114    procedure getrf
 115      (M     : Natural;
 116       N     : Natural;
 117       A     : in out Real_Matrix;
 118       Ld_A  : Positive;
 119       I_Piv : out Integer_Vector;
 120       Info  : access Integer)
 121    is
 122    begin
 123       if Is_Real then
 124          declare
 125             type A_Ptr is
 126                access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
 127             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 128          begin
 129             sgetrf (M, N, Conv_A (A'Address).all, Ld_A,
 130                     LAPACK.Integer_Vector (I_Piv), Info);
 131          end;
 132
 133       elsif Is_Double_Precision then
 134          declare
 135             type A_Ptr is
 136                access all Double_Precision_Matrix (A'Range (1), A'Range (2));
 137             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 138          begin
 139             dgetrf (M, N, Conv_A (A'Address).all, Ld_A,
 140                     LAPACK.Integer_Vector (I_Piv), Info);
 141          end;
 142
 143       else
 144          declare
 145             DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
 146          begin
 147             DP_A := To_Double_Precision (A);
 148             dgetrf (M, N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), Info);
 149             A := To_Real (DP_A);
 150          end;
 151       end if;
 152    end getrf;
 153
 154    -----------
 155    -- getri --
 156    -----------
 157
 158    procedure getri
 159      (N      : Natural;
 160       A      : in out Real_Matrix;
 161       Ld_A   : Positive;
 162       I_Piv  : Integer_Vector;
 163       Work   : in out Real_Vector;
 164       L_Work : Integer;
 165       Info   : access Integer)
 166    is
 167    begin
 168       if Is_Real then
 169          declare
 170             type A_Ptr is
 171                access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
 172             type Work_Ptr is
 173                access all BLAS.Real_Vector (Work'Range);
 174             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 175             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 176          begin
 177             sgetri (N, Conv_A (A'Address).all, Ld_A,
 178                     LAPACK.Integer_Vector (I_Piv),
 179                     Conv_Work (Work'Address).all, L_Work,
 180                     Info);
 181          end;
 182
 183       elsif Is_Double_Precision then
 184          declare
 185             type A_Ptr is
 186                access all Double_Precision_Matrix (A'Range (1), A'Range (2));
 187             type Work_Ptr is
 188                access all Double_Precision_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             dgetri (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       else
 199          declare
 200             DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
 201             DP_Work : Double_Precision_Vector (Work'Range);
 202          begin
 203             DP_A := To_Double_Precision (A);
 204             dgetri (N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv),
 205                     DP_Work, L_Work, Info);
 206             A := To_Real (DP_A);
 207             Work (1) := To_Real (DP_Work (1));
 208          end;
 209       end if;
 210    end getri;
 211
 212    -----------
 213    -- getrs --
 214    -----------
 215
 216    procedure getrs
 217      (Trans  : access constant Character;
 218       N      : Natural;
 219       N_Rhs  : Natural;
 220       A      : Real_Matrix;
 221       Ld_A   : Positive;
 222       I_Piv  : Integer_Vector;
 223       B      : in out Real_Matrix;
 224       Ld_B   : Positive;
 225       Info   : access Integer)
 226    is
 227    begin
 228       if Is_Real then
 229          declare
 230             subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2));
 231             type B_Ptr is
 232                access all BLAS.Real_Matrix (B'Range (1), B'Range (2));
 233             function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
 234             function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
 235          begin
 236             sgetrs (Trans, N, N_Rhs,
 237                     Conv_A (A), Ld_A,
 238                     LAPACK.Integer_Vector (I_Piv),
 239                     Conv_B (B'Address).all, Ld_B,
 240                     Info);
 241          end;
 242
 243       elsif Is_Double_Precision then
 244          declare
 245             subtype A_Type is
 246                Double_Precision_Matrix (A'Range (1), A'Range (2));
 247             type B_Ptr is
 248                access all Double_Precision_Matrix (B'Range (1), B'Range (2));
 249             function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type);
 250             function Conv_B is new Unchecked_Conversion (Address, B_Ptr);
 251          begin
 252             dgetrs (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       else
 260          declare
 261             DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
 262             DP_B : Double_Precision_Matrix (B'Range (1), B'Range (2));
 263          begin
 264             DP_A := To_Double_Precision (A);
 265             DP_B := To_Double_Precision (B);
 266             dgetrs (Trans, N, N_Rhs,
 267                     DP_A, Ld_A,
 268                     LAPACK.Integer_Vector (I_Piv),
 269                     DP_B, Ld_B,
 270                     Info);
 271             B := To_Real (DP_B);
 272          end;
 273       end if;
 274    end getrs;
 275
 276    -----------
 277    -- orgtr --
 278    -----------
 279
 280    procedure orgtr
 281      (Uplo   : access constant Character;
 282       N      : Natural;
 283       A      : in out Real_Matrix;
 284       Ld_A   : Positive;
 285       Tau    : in Real_Vector;
 286       Work   : out Real_Vector;
 287       L_Work : Integer;
 288       Info   : access Integer)
 289    is
 290    begin
 291       if Is_Real then
 292          declare
 293             type A_Ptr is
 294                access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
 295             subtype Tau_Type is BLAS.Real_Vector (Tau'Range);
 296             type Work_Ptr is
 297                access all BLAS.Real_Vector (Work'Range);
 298             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 299             function Conv_Tau is
 300                new Unchecked_Conversion (Real_Vector, Tau_Type);
 301             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 302          begin
 303             sorgtr (Uplo, N,
 304                     Conv_A (A'Address).all, Ld_A,
 305                     Conv_Tau (Tau),
 306                     Conv_Work (Work'Address).all, L_Work,
 307                     Info);
 308          end;
 309
 310       elsif Is_Double_Precision then
 311          declare
 312             type A_Ptr is
 313                access all Double_Precision_Matrix (A'Range (1), A'Range (2));
 314             subtype Tau_Type is Double_Precision_Vector (Tau'Range);
 315             type Work_Ptr is
 316                access all Double_Precision_Vector (Work'Range);
 317             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 318             function Conv_Tau is
 319                new Unchecked_Conversion (Real_Vector, Tau_Type);
 320             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 321          begin
 322             dorgtr (Uplo, N,
 323                     Conv_A (A'Address).all, Ld_A,
 324                     Conv_Tau (Tau),
 325                     Conv_Work (Work'Address).all, L_Work,
 326                     Info);
 327          end;
 328
 329       else
 330          declare
 331             DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2));
 332             DP_Work : Double_Precision_Vector (Work'Range);
 333             DP_Tau  : Double_Precision_Vector (Tau'Range);
 334          begin
 335             DP_A := To_Double_Precision (A);
 336             DP_Tau := To_Double_Precision (Tau);
 337             dorgtr (Uplo, N, DP_A, Ld_A, DP_Tau, DP_Work, L_Work, Info);
 338             A := To_Real (DP_A);
 339             Work (1) := To_Real (DP_Work (1));
 340          end;
 341       end if;
 342    end orgtr;
 343
 344    -----------
 345    -- steqr --
 346    -----------
 347
 348    procedure steqr
 349      (Comp_Z : access constant Character;
 350       N      : Natural;
 351       D      : in out Real_Vector;
 352       E      : in out Real_Vector;
 353       Z      : in out Real_Matrix;
 354       Ld_Z   : Positive;
 355       Work   : out Real_Vector;
 356       Info   : access Integer)
 357    is
 358    begin
 359       if Is_Real then
 360          declare
 361             type D_Ptr is access all BLAS.Real_Vector (D'Range);
 362             type E_Ptr is access all BLAS.Real_Vector (E'Range);
 363             type Z_Ptr is
 364                access all BLAS.Real_Matrix (Z'Range (1), Z'Range (2));
 365             type Work_Ptr is
 366                access all BLAS.Real_Vector (Work'Range);
 367             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 368             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 369             function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
 370             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 371          begin
 372             ssteqr (Comp_Z, N,
 373                     Conv_D (D'Address).all,
 374                     Conv_E (E'Address).all,
 375                     Conv_Z (Z'Address).all,
 376                     Ld_Z,
 377                     Conv_Work (Work'Address).all,
 378                     Info);
 379          end;
 380
 381       elsif Is_Double_Precision then
 382          declare
 383             type D_Ptr is access all Double_Precision_Vector (D'Range);
 384             type E_Ptr is access all Double_Precision_Vector (E'Range);
 385             type Z_Ptr is
 386                access all Double_Precision_Matrix (Z'Range (1), Z'Range (2));
 387             type Work_Ptr is
 388                access all Double_Precision_Vector (Work'Range);
 389             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 390             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 391             function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr);
 392             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 393          begin
 394             dsteqr (Comp_Z, N,
 395                     Conv_D (D'Address).all,
 396                     Conv_E (E'Address).all,
 397                     Conv_Z (Z'Address).all,
 398                     Ld_Z,
 399                     Conv_Work (Work'Address).all,
 400                     Info);
 401          end;
 402
 403       else
 404          declare
 405             DP_D    : Double_Precision_Vector (D'Range);
 406             DP_E    : Double_Precision_Vector (E'Range);
 407             DP_Z    : Double_Precision_Matrix (Z'Range (1), Z'Range (2));
 408             DP_Work : Double_Precision_Vector (Work'Range);
 409          begin
 410             DP_D := To_Double_Precision (D);
 411             DP_E := To_Double_Precision (E);
 412
 413             if Comp_Z.all = 'V' then
 414                DP_Z := To_Double_Precision (Z);
 415             end if;
 416
 417             dsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info);
 418
 419             D := To_Real (DP_D);
 420             E := To_Real (DP_E);
 421             Z := To_Real (DP_Z);
 422          end;
 423       end if;
 424    end steqr;
 425
 426    -----------
 427    -- sterf --
 428    -----------
 429
 430    procedure sterf
 431      (N    : Natural;
 432       D    : in out Real_Vector;
 433       E    : in out Real_Vector;
 434       Info : access Integer)
 435    is
 436    begin
 437       if Is_Real then
 438          declare
 439             type D_Ptr is access all BLAS.Real_Vector (D'Range);
 440             type E_Ptr is access all BLAS.Real_Vector (E'Range);
 441             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 442             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 443          begin
 444             ssterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info);
 445          end;
 446
 447       elsif Is_Double_Precision 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             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 452             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 453          begin
 454             dsterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info);
 455          end;
 456
 457       else
 458          declare
 459             DP_D    : Double_Precision_Vector (D'Range);
 460             DP_E    : Double_Precision_Vector (E'Range);
 461
 462          begin
 463             DP_D := To_Double_Precision (D);
 464             DP_E := To_Double_Precision (E);
 465
 466             dsterf (N, DP_D, DP_E, Info);
 467
 468             D := To_Real (DP_D);
 469             E := To_Real (DP_E);
 470          end;
 471       end if;
 472    end sterf;
 473
 474    -----------
 475    -- sytrd --
 476    -----------
 477
 478    procedure sytrd
 479      (Uplo   : access constant Character;
 480       N      : Natural;
 481       A      : in out Real_Matrix;
 482       Ld_A   : Positive;
 483       D      : out Real_Vector;
 484       E      : out Real_Vector;
 485       Tau    : out Real_Vector;
 486       Work   : out Real_Vector;
 487       L_Work : Integer;
 488       Info   : access Integer)
 489    is
 490    begin
 491       if Is_Real then
 492          declare
 493             type A_Ptr is
 494                access all BLAS.Real_Matrix (A'Range (1), A'Range (2));
 495             type D_Ptr is access all BLAS.Real_Vector (D'Range);
 496             type E_Ptr is access all BLAS.Real_Vector (E'Range);
 497             type Tau_Ptr is access all BLAS.Real_Vector (Tau'Range);
 498             type Work_Ptr is
 499                access all BLAS.Real_Vector (Work'Range);
 500             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 501             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 502             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 503             function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr);
 504             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 505          begin
 506             ssytrd (Uplo, N,
 507                     Conv_A (A'Address).all, Ld_A,
 508                     Conv_D (D'Address).all,
 509                     Conv_E (E'Address).all,
 510                     Conv_Tau (Tau'Address).all,
 511                     Conv_Work (Work'Address).all,
 512                     L_Work,
 513                     Info);
 514          end;
 515
 516       elsif Is_Double_Precision then
 517          declare
 518             type A_Ptr is
 519                access all Double_Precision_Matrix (A'Range (1), A'Range (2));
 520             type D_Ptr is access all Double_Precision_Vector (D'Range);
 521             type E_Ptr is access all Double_Precision_Vector (E'Range);
 522             type Tau_Ptr is access all Double_Precision_Vector (Tau'Range);
 523             type Work_Ptr is
 524                access all Double_Precision_Vector (Work'Range);
 525             function Conv_A is new Unchecked_Conversion (Address, A_Ptr);
 526             function Conv_D is new Unchecked_Conversion (Address, D_Ptr);
 527             function Conv_E is new Unchecked_Conversion (Address, E_Ptr);
 528             function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr);
 529             function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr);
 530          begin
 531             dsytrd (Uplo, N,
 532                     Conv_A (A'Address).all, Ld_A,
 533                     Conv_D (D'Address).all,
 534                     Conv_E (E'Address).all,
 535                     Conv_Tau (Tau'Address).all,
 536                     Conv_Work (Work'Address).all,
 537                     L_Work,
 538                     Info);
 539          end;
 540
 541       else
 542          declare
 543             DP_A    : Double_Precision_Matrix (A'Range (1), A'Range (2));
 544             DP_D    : Double_Precision_Vector (D'Range);
 545             DP_E    : Double_Precision_Vector (E'Range);
 546             DP_Tau  : Double_Precision_Vector (Tau'Range);
 547             DP_Work : Double_Precision_Vector (Work'Range);
 548          begin
 549             DP_A := To_Double_Precision (A);
 550
 551             dsytrd (Uplo, N, DP_A, Ld_A, DP_D, DP_E, DP_Tau,
 552                     DP_Work, L_Work, Info);
 553
 554             if L_Work /= -1 then
 555                A := To_Real (DP_A);
 556                D := To_Real (DP_D);
 557                E := To_Real (DP_E);
 558                Tau := To_Real (DP_Tau);
 559             end if;
 560
 561             Work (1) := To_Real (DP_Work (1));
 562          end;
 563       end if;
 564    end sytrd;
 565
 566 end System.Generic_Real_LAPACK;