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