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