gcc/ada/a-ngrear.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                     ADA.NUMERICS.GENERIC_REAL_ARRAYS                     --
   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 System; use System;
  35 with System.Generic_Real_BLAS;
  36 with System.Generic_Real_LAPACK;
  37 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
  38
  39 package body Ada.Numerics.Generic_Real_Arrays is
  40
  41    --  Operations involving inner products use BLAS library implementations.
  42    --  This allows larger matrices and vectors to be computed efficiently,
  43    --  taking into account memory hierarchy issues and vector instructions
  44    --  that vary widely between machines.
  45
  46    --  Operations that are defined in terms of operations on the type Real,
  47    --  such as addition, subtraction and scaling, are computed in the canonical
  48    --  way looping over all elements.
  49
  50    --  Operations for solving linear systems and computing determinant,
  51    --  eigenvalues, eigensystem and inverse, are implemented using the
  52    --  LAPACK library.
  53
  54    package BLAS is
  55       new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);
  56
  57    package LAPACK is
  58       new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);
  59
  60    use BLAS, LAPACK;
  61
  62    --  Procedure versions of functions returning unconstrained values.
  63    --  This allows for inlining the function wrapper.
  64
  65    procedure Eigenvalues (A : Real_Matrix; Values : out Real_Vector);
  66    procedure Inverse   (A : Real_Matrix; R : out Real_Matrix);
  67    procedure Solve     (A : Real_Matrix; X : Real_Vector; B : out Real_Vector);
  68    procedure Solve     (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix);
  69
  70    procedure Transpose is new
  71      Generic_Array_Operations.Transpose
  72        (Scalar        => Real'Base,
  73         Matrix        => Real_Matrix);
  74
  75    --  Helper function that raises a Constraint_Error is the argument is
  76    --  not a square matrix, and otherwise returns its length.
  77
  78    function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
  79
  80    --  Instantiating the following subprograms directly would lead to
  81    --  name clashes, so use a local package.
  82
  83    package Instantiations is
  84
  85       function "+" is new
  86         Vector_Elementwise_Operation
  87           (X_Scalar      => Real'Base,
  88            Result_Scalar => Real'Base,
  89            X_Vector      => Real_Vector,
  90            Result_Vector => Real_Vector,
  91            Operation     => "+");
  92
  93       function "+" is new
  94         Matrix_Elementwise_Operation
  95           (X_Scalar      => Real'Base,
  96            Result_Scalar => Real'Base,
  97            X_Matrix      => Real_Matrix,
  98            Result_Matrix => Real_Matrix,
  99            Operation     => "+");
 100
 101       function "+" is new
 102         Vector_Vector_Elementwise_Operation
 103           (Left_Scalar   => Real'Base,
 104            Right_Scalar  => Real'Base,
 105            Result_Scalar => Real'Base,
 106            Left_Vector   => Real_Vector,
 107            Right_Vector  => Real_Vector,
 108            Result_Vector => Real_Vector,
 109            Operation     => "+");
 110
 111       function "+" is new
 112         Matrix_Matrix_Elementwise_Operation
 113           (Left_Scalar   => Real'Base,
 114            Right_Scalar  => Real'Base,
 115            Result_Scalar => Real'Base,
 116            Left_Matrix   => Real_Matrix,
 117            Right_Matrix  => Real_Matrix,
 118            Result_Matrix => Real_Matrix,
 119            Operation     => "+");
 120
 121       function "-" is new
 122         Vector_Elementwise_Operation
 123           (X_Scalar      => Real'Base,
 124            Result_Scalar => Real'Base,
 125            X_Vector      => Real_Vector,
 126            Result_Vector => Real_Vector,
 127            Operation     => "-");
 128
 129       function "-" is new
 130         Matrix_Elementwise_Operation
 131           (X_Scalar      => Real'Base,
 132            Result_Scalar => Real'Base,
 133            X_Matrix      => Real_Matrix,
 134            Result_Matrix => Real_Matrix,
 135            Operation     => "-");
 136
 137       function "-" is new
 138         Vector_Vector_Elementwise_Operation
 139           (Left_Scalar   => Real'Base,
 140            Right_Scalar  => Real'Base,
 141            Result_Scalar => Real'Base,
 142            Left_Vector   => Real_Vector,
 143            Right_Vector  => Real_Vector,
 144            Result_Vector => Real_Vector,
 145            Operation     => "-");
 146
 147       function "-" is new
 148         Matrix_Matrix_Elementwise_Operation
 149           (Left_Scalar   => Real'Base,
 150            Right_Scalar  => Real'Base,
 151            Result_Scalar => Real'Base,
 152            Left_Matrix   => Real_Matrix,
 153            Right_Matrix  => Real_Matrix,
 154            Result_Matrix => Real_Matrix,
 155            Operation     => "-");
 156
 157       function "*" is new
 158         Scalar_Vector_Elementwise_Operation
 159           (Left_Scalar   => Real'Base,
 160            Right_Scalar  => Real'Base,
 161            Result_Scalar => Real'Base,
 162            Right_Vector  => Real_Vector,
 163            Result_Vector => Real_Vector,
 164            Operation     => "*");
 165
 166       function "*" is new
 167         Scalar_Matrix_Elementwise_Operation
 168           (Left_Scalar   => Real'Base,
 169            Right_Scalar  => Real'Base,
 170            Result_Scalar => Real'Base,
 171            Right_Matrix  => Real_Matrix,
 172            Result_Matrix => Real_Matrix,
 173            Operation     => "*");
 174
 175       function "*" is new
 176         Vector_Scalar_Elementwise_Operation
 177           (Left_Scalar   => Real'Base,
 178            Right_Scalar  => Real'Base,
 179            Result_Scalar => Real'Base,
 180            Left_Vector   => Real_Vector,
 181            Result_Vector => Real_Vector,
 182            Operation     => "*");
 183
 184       function "*" is new
 185         Matrix_Scalar_Elementwise_Operation
 186           (Left_Scalar   => Real'Base,
 187            Right_Scalar  => Real'Base,
 188            Result_Scalar => Real'Base,
 189            Left_Matrix   => Real_Matrix,
 190            Result_Matrix => Real_Matrix,
 191            Operation     => "*");
 192
 193       function "*" is new
 194         Outer_Product
 195           (Left_Scalar   => Real'Base,
 196            Right_Scalar  => Real'Base,
 197            Result_Scalar => Real'Base,
 198            Left_Vector   => Real_Vector,
 199            Right_Vector  => Real_Vector,
 200            Matrix        => Real_Matrix);
 201
 202       function "/" is new
 203         Vector_Scalar_Elementwise_Operation
 204           (Left_Scalar   => Real'Base,
 205            Right_Scalar  => Real'Base,
 206            Result_Scalar => Real'Base,
 207            Left_Vector   => Real_Vector,
 208            Result_Vector => Real_Vector,
 209            Operation     => "/");
 210
 211       function "/" is new
 212         Matrix_Scalar_Elementwise_Operation
 213           (Left_Scalar   => Real'Base,
 214            Right_Scalar  => Real'Base,
 215            Result_Scalar => Real'Base,
 216            Left_Matrix   => Real_Matrix,
 217            Result_Matrix => Real_Matrix,
 218            Operation     => "/");
 219
 220       function "abs" is new
 221         Vector_Elementwise_Operation
 222           (X_Scalar      => Real'Base,
 223            Result_Scalar => Real'Base,
 224            X_Vector      => Real_Vector,
 225            Result_Vector => Real_Vector,
 226            Operation     => "abs");
 227
 228       function "abs" is new
 229         Matrix_Elementwise_Operation
 230           (X_Scalar      => Real'Base,
 231            Result_Scalar => Real'Base,
 232            X_Matrix      => Real_Matrix,
 233            Result_Matrix => Real_Matrix,
 234            Operation     => "abs");
 235
 236       function Unit_Matrix is new
 237         Generic_Array_Operations.Unit_Matrix
 238           (Scalar        => Real'Base,
 239            Matrix        => Real_Matrix,
 240            Zero          => 0.0,
 241            One           => 1.0);
 242
 243       function Unit_Vector is new
 244         Generic_Array_Operations.Unit_Vector
 245           (Scalar        => Real'Base,
 246            Vector        => Real_Vector,
 247            Zero          => 0.0,
 248            One           => 1.0);
 249
 250    end Instantiations;
 251
 252    ---------
 253    -- "+" --
 254    ---------
 255
 256    function "+" (Right : Real_Vector) return Real_Vector
 257       renames Instantiations."+";
 258
 259    function "+" (Right : Real_Matrix) return Real_Matrix
 260       renames Instantiations."+";
 261
 262    function "+" (Left, Right : Real_Vector) return Real_Vector
 263       renames Instantiations."+";
 264
 265    function "+" (Left, Right : Real_Matrix) return Real_Matrix
 266       renames Instantiations."+";
 267
 268    ---------
 269    -- "-" --
 270    ---------
 271
 272    function "-" (Right : Real_Vector) return Real_Vector
 273       renames Instantiations."-";
 274
 275    function "-" (Right : Real_Matrix) return Real_Matrix
 276       renames Instantiations."-";
 277
 278    function "-" (Left, Right : Real_Vector) return Real_Vector
 279       renames Instantiations."-";
 280
 281    function "-" (Left, Right : Real_Matrix) return Real_Matrix
 282       renames Instantiations."-";
 283
 284    ---------
 285    -- "*" --
 286    ---------
 287
 288    --  Scalar multiplication
 289
 290    function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector
 291       renames Instantiations."*";
 292
 293    function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector
 294       renames Instantiations."*";
 295
 296    function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix
 297       renames Instantiations."*";
 298
 299    function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
 300       renames Instantiations."*";
 301
 302    --  Vector multiplication
 303
 304    function "*" (Left, Right : Real_Vector) return Real'Base is
 305    begin
 306       if Left'Length /= Right'Length then
 307          raise Constraint_Error with
 308             "vectors are of different length in inner product";
 309       end if;
 310
 311       return dot (Left'Length, X => Left, Y => Right);
 312    end "*";
 313
 314    function "*" (Left, Right : Real_Vector) return Real_Matrix
 315       renames Instantiations."*";
 316
 317    function "*"
 318      (Left : Real_Vector;
 319       Right : Real_Matrix) return Real_Vector
 320    is
 321       R : Real_Vector (Right'Range (2));
 322
 323    begin
 324       if Left'Length /= Right'Length (1) then
 325          raise Constraint_Error with
 326            "incompatible dimensions in vector-matrix multiplication";
 327       end if;
 328
 329       gemv (Trans => No_Trans'Access,
 330             M     => Right'Length (2),
 331             N     => Right'Length (1),
 332             A     => Right,
 333             Ld_A  => Right'Length (2),
 334             X     => Left,
 335             Y     => R);
 336
 337       return R;
 338    end "*";
 339
 340    function "*"
 341      (Left : Real_Matrix;
 342       Right : Real_Vector) return Real_Vector
 343    is
 344       R : Real_Vector (Left'Range (1));
 345
 346    begin
 347       if Left'Length (2) /= Right'Length then
 348          raise Constraint_Error with
 349             "incompatible dimensions in matrix-vector multiplication";
 350       end if;
 351
 352       gemv (Trans => Trans'Access,
 353             M     => Left'Length (2),
 354             N     => Left'Length (1),
 355             A     => Left,
 356             Ld_A  => Left'Length (2),
 357             X     => Right,
 358             Y     => R);
 359
 360       return R;
 361    end "*";
 362
 363    --  Matrix Multiplication
 364
 365    function "*" (Left, Right : Real_Matrix) return Real_Matrix is
 366       R : Real_Matrix (Left'Range (1), Right'Range (2));
 367
 368    begin
 369       if Left'Length (2) /= Right'Length (1) then
 370          raise Constraint_Error with
 371             "incompatible dimensions in matrix-matrix multipication";
 372       end if;
 373
 374       gemm (Trans_A => No_Trans'Access,
 375             Trans_B => No_Trans'Access,
 376             M       => Right'Length (2),
 377             N       => Left'Length (1),
 378             K       => Right'Length (1),
 379             A       => Right,
 380             Ld_A    => Right'Length (2),
 381             B       => Left,
 382             Ld_B    => Left'Length (2),
 383             C       => R,
 384             Ld_C    => R'Length (2));
 385
 386       return R;
 387    end "*";
 388
 389    ---------
 390    -- "/" --
 391    ---------
 392
 393    function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector
 394       renames Instantiations."/";
 395
 396    function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix
 397       renames Instantiations."/";
 398
 399    -----------
 400    -- "abs" --
 401    -----------
 402
 403    function "abs" (Right : Real_Vector) return Real'Base is
 404    begin
 405       return nrm2 (Right'Length, Right);
 406    end "abs";
 407
 408    function "abs" (Right : Real_Vector) return Real_Vector
 409       renames Instantiations."abs";
 410
 411    function "abs" (Right : Real_Matrix) return Real_Matrix
 412       renames Instantiations."abs";
 413
 414    -----------------
 415    -- Determinant --
 416    -----------------
 417
 418    function Determinant (A : Real_Matrix) return Real'Base is
 419       N    : constant Integer := Length (A);
 420       LU   : Real_Matrix (1 .. N, 1 .. N) := A;
 421       Piv  : Integer_Vector (1 .. N);
 422       Info : aliased Integer := -1;
 423       Det  : Real := 1.0;
 424
 425    begin
 426       getrf (M     => N,
 427              N     => N,
 428              A     => LU,
 429              Ld_A  => N,
 430              I_Piv => Piv,
 431              Info  => Info'Access);
 432
 433       if Info /= 0 then
 434          raise Constraint_Error with "ill-conditioned matrix";
 435       end if;
 436
 437       for J in 1 .. N loop
 438          if Piv (J) /= J then
 439             Det := -Det * LU (J, J);
 440          else
 441             Det := Det * LU (J, J);
 442          end if;
 443       end loop;
 444
 445       return Det;
 446    end Determinant;
 447
 448    -----------------
 449    -- Eigensystem --
 450    -----------------
 451
 452    procedure Eigensystem
 453      (A       : Real_Matrix;
 454       Values  : out Real_Vector;
 455       Vectors : out Real_Matrix)
 456    is
 457       N      : constant Natural := Length (A);
 458       E      : Real_Vector (1 .. N);
 459       Tau    : Real_Vector (1 .. N);
 460       L_Work : Real_Vector (1 .. 1);
 461       Info   : aliased Integer;
 462
 463    begin
 464       if Values'Length /= N then
 465          raise Constraint_Error with "wrong length for output vector";
 466       end if;
 467
 468       if N = 0 then
 469          return;
 470       end if;
 471
 472       --  Initialize working matrix and check for symmetric input matrix
 473
 474       Transpose (A, Vectors);
 475
 476       if A /= Vectors then
 477          raise Argument_Error with "matrix not symmetric";
 478       end if;
 479
 480       --  Compute size of additional working space
 481
 482       sytrd (Uplo   => Lower'Access,
 483              N      => N,
 484              A      => Vectors,
 485              Ld_A   => N,
 486              D      => Values,
 487              E      => E,
 488              Tau    => Tau,
 489              Work   => L_Work,
 490              L_Work => -1,
 491              Info   => Info'Access);
 492
 493       declare
 494          Work   : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
 495          Comp_Z : aliased constant Character := 'V';
 496
 497       begin
 498          --  Reduce matrix to tridiagonal form
 499
 500          sytrd (Uplo   => Lower'Access,
 501                 N      => N,
 502                 A      => Vectors,
 503                 Ld_A   => A'Length (1),
 504                 D      => Values,
 505                 E      => E,
 506                 Tau    => Tau,
 507                 Work   => Work,
 508                 L_Work => Work'Length,
 509                 Info   => Info'Access);
 510
 511          if Info /= 0 then
 512             raise Program_Error;
 513          end if;
 514
 515          --  Generate the real orthogonal matrix determined by sytrd
 516
 517          orgtr (Uplo   => Lower'Access,
 518                 N      => N,
 519                 A      => Vectors,
 520                 Ld_A   => N,
 521                 Tau    => Tau,
 522                 Work   => Work,
 523                 L_Work => Work'Length,
 524                 Info   => Info'Access);
 525
 526          if Info /= 0 then
 527             raise Program_Error;
 528          end if;
 529
 530          --  Compute all eigenvalues and eigenvectors using QR algorithm
 531
 532          steqr (Comp_Z => Comp_Z'Access,
 533                 N      => N,
 534                 D      => Values,
 535                 E      => E,
 536                 Z      => Vectors,
 537                 Ld_Z   => N,
 538                 Work   => Work,
 539                 Info   => Info'Access);
 540
 541          if Info /= 0 then
 542             raise Constraint_Error with
 543                "eigensystem computation failed to converge";
 544          end if;
 545       end;
 546    end Eigensystem;
 547
 548    -----------------
 549    -- Eigenvalues --
 550    -----------------
 551
 552    procedure Eigenvalues
 553      (A      : Real_Matrix;
 554       Values : out Real_Vector)
 555    is
 556       N      : constant Natural := Length (A);
 557       B      : Real_Matrix (1 .. N, 1 .. N);
 558       E      : Real_Vector (1 .. N);
 559       Tau    : Real_Vector (1 .. N);
 560       L_Work : Real_Vector (1 .. 1);
 561       Info   : aliased Integer;
 562
 563    begin
 564       if Values'Length /= N then
 565          raise Constraint_Error with "wrong length for output vector";
 566       end if;
 567
 568       if N = 0 then
 569          return;
 570       end if;
 571
 572       --  Initialize working matrix and check for symmetric input matrix
 573
 574       Transpose (A, B);
 575
 576       if A /= B then
 577          raise Argument_Error with "matrix not symmetric";
 578       end if;
 579
 580       --  Find size of work area
 581
 582       sytrd (Uplo   => Lower'Access,
 583              N      => N,
 584              A      => B,
 585              Ld_A   => N,
 586              D      => Values,
 587              E      => E,
 588              Tau    => Tau,
 589              Work   => L_Work,
 590              L_Work => -1,
 591              Info   => Info'Access);
 592
 593       declare
 594          Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
 595
 596       begin
 597          --  Reduce matrix to tridiagonal form
 598
 599          sytrd (Uplo   => Lower'Access,
 600                 N      => N,
 601                 A      => B,
 602                 Ld_A   => A'Length (1),
 603                 D      => Values,
 604                 E      => E,
 605                 Tau    => Tau,
 606                 Work   => Work,
 607                 L_Work => Work'Length,
 608                 Info   => Info'Access);
 609
 610          if Info /= 0 then
 611             raise Constraint_Error;
 612          end if;
 613
 614          --  Compute all eigenvalues using QR algorithm
 615
 616          sterf (N      => N,
 617                 D      => Values,
 618                 E      => E,
 619                 Info   => Info'Access);
 620
 621          if Info /= 0 then
 622             raise Constraint_Error with
 623                "eigenvalues computation failed to converge";
 624          end if;
 625       end;
 626    end Eigenvalues;
 627
 628    function Eigenvalues (A : Real_Matrix) return Real_Vector is
 629       R : Real_Vector (A'Range (1));
 630    begin
 631       Eigenvalues (A, R);
 632       return R;
 633    end Eigenvalues;
 634
 635    -------------
 636    -- Inverse --
 637    -------------
 638
 639    procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is
 640       N      : constant Integer := Length (A);
 641       Piv    : Integer_Vector (1 .. N);
 642       L_Work : Real_Vector (1 .. 1);
 643       Info   : aliased Integer := -1;
 644
 645    begin
 646       --  All computations are done using column-major order, but this works
 647       --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
 648
 649       R := A;
 650
 651       --  Compute LU decomposition
 652
 653       getrf (M      => N,
 654              N      => N,
 655              A      => R,
 656              Ld_A   => N,
 657              I_Piv  => Piv,
 658              Info   => Info'Access);
 659
 660       if Info /= 0 then
 661          raise Constraint_Error with "inverting singular matrix";
 662       end if;
 663
 664       --  Determine size of work area
 665
 666       getri (N      => N,
 667              A      => R,
 668              Ld_A   => N,
 669              I_Piv  => Piv,
 670              Work   => L_Work,
 671              L_Work => -1,
 672              Info   => Info'Access);
 673
 674       if Info /= 0 then
 675          raise Constraint_Error;
 676       end if;
 677
 678       declare
 679          Work : Real_Vector (1 .. Integer (L_Work (1)));
 680       begin
 681          --  Compute inverse from LU decomposition
 682
 683          getri (N      => N,
 684                 A      => R,
 685                 Ld_A   => N,
 686                 I_Piv  => Piv,
 687                 Work   => Work,
 688                 L_Work => Work'Length,
 689                 Info   => Info'Access);
 690
 691          if Info /= 0 then
 692             raise Constraint_Error with "inverting singular matrix";
 693          end if;
 694
 695          --  ??? Should iterate with gerfs, based on implementation advice
 696       end;
 697    end Inverse;
 698
 699    function Inverse (A : Real_Matrix) return Real_Matrix is
 700       R : Real_Matrix (A'Range (2), A'Range (1));
 701    begin
 702       Inverse (A, R);
 703       return R;
 704    end Inverse;
 705
 706    -----------
 707    -- Solve --
 708    -----------
 709
 710    procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector) is
 711    begin
 712       if Length (A) /= X'Length then
 713          raise Constraint_Error with
 714            "incompatible matrix and vector dimensions";
 715       end if;
 716
 717       --  ??? Should solve directly, is faster and more accurate
 718
 719       B := Inverse (A) * X;
 720    end Solve;
 721
 722    procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix) is
 723    begin
 724       if Length (A) /= X'Length (1) then
 725          raise Constraint_Error with "incompatible matrix dimensions";
 726       end if;
 727
 728       --  ??? Should solve directly, is faster and more accurate
 729
 730       B := Inverse (A) * X;
 731    end Solve;
 732
 733    function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
 734       B : Real_Vector (A'Range (2));
 735    begin
 736       Solve (A, X, B);
 737       return B;
 738    end Solve;
 739
 740    function Solve (A, X : Real_Matrix) return Real_Matrix is
 741       B : Real_Matrix (A'Range (2), X'Range (2));
 742    begin
 743       Solve (A, X, B);
 744       return B;
 745    end Solve;
 746
 747    ---------------
 748    -- Transpose --
 749    ---------------
 750
 751    function Transpose (X : Real_Matrix) return Real_Matrix is
 752       R : Real_Matrix (X'Range (2), X'Range (1));
 753    begin
 754       Transpose (X, R);
 755
 756       return R;
 757    end Transpose;
 758
 759    -----------------
 760    -- Unit_Matrix --
 761    -----------------
 762
 763    function Unit_Matrix
 764      (Order   : Positive;
 765       First_1 : Integer := 1;
 766       First_2 : Integer := 1) return Real_Matrix
 767      renames Instantiations.Unit_Matrix;
 768
 769    -----------------
 770    -- Unit_Vector --
 771    -----------------
 772
 773    function Unit_Vector
 774      (Index : Integer;
 775       Order : Positive;
 776       First : Integer := 1) return Real_Vector
 777      renames Instantiations.Unit_Vector;
 778
 779 end Ada.Numerics.Generic_Real_Arrays;