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-2007, 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       Tau    : Real_Vector (1 .. N);
 459       L_Work : Real_Vector (1 .. 1);
 460       Info   : aliased Integer;
 461
 462       E : Real_Vector (1 .. N);
 463       pragma Warnings (Off, E);
 464
 465    begin
 466       if Values'Length /= N then
 467          raise Constraint_Error with "wrong length for output vector";
 468       end if;
 469
 470       if N = 0 then
 471          return;
 472       end if;
 473
 474       --  Initialize working matrix and check for symmetric input matrix
 475
 476       Transpose (A, Vectors);
 477
 478       if A /= Vectors then
 479          raise Argument_Error with "matrix not symmetric";
 480       end if;
 481
 482       --  Compute size of additional working space
 483
 484       sytrd (Uplo   => Lower'Access,
 485              N      => N,
 486              A      => Vectors,
 487              Ld_A   => N,
 488              D      => Values,
 489              E      => E,
 490              Tau    => Tau,
 491              Work   => L_Work,
 492              L_Work => -1,
 493              Info   => Info'Access);
 494
 495       declare
 496          Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
 497          pragma Warnings (Off, Work);
 498
 499          Comp_Z : aliased constant Character := 'V';
 500
 501       begin
 502          --  Reduce matrix to tridiagonal form
 503
 504          sytrd (Uplo   => Lower'Access,
 505                 N      => N,
 506                 A      => Vectors,
 507                 Ld_A   => A'Length (1),
 508                 D      => Values,
 509                 E      => E,
 510                 Tau    => Tau,
 511                 Work   => Work,
 512                 L_Work => Work'Length,
 513                 Info   => Info'Access);
 514
 515          if Info /= 0 then
 516             raise Program_Error;
 517          end if;
 518
 519          --  Generate the real orthogonal matrix determined by sytrd
 520
 521          orgtr (Uplo   => Lower'Access,
 522                 N      => N,
 523                 A      => Vectors,
 524                 Ld_A   => N,
 525                 Tau    => Tau,
 526                 Work   => Work,
 527                 L_Work => Work'Length,
 528                 Info   => Info'Access);
 529
 530          if Info /= 0 then
 531             raise Program_Error;
 532          end if;
 533
 534          --  Compute all eigenvalues and eigenvectors using QR algorithm
 535
 536          steqr (Comp_Z => Comp_Z'Access,
 537                 N      => N,
 538                 D      => Values,
 539                 E      => E,
 540                 Z      => Vectors,
 541                 Ld_Z   => N,
 542                 Work   => Work,
 543                 Info   => Info'Access);
 544
 545          if Info /= 0 then
 546             raise Constraint_Error with
 547                "eigensystem computation failed to converge";
 548          end if;
 549       end;
 550    end Eigensystem;
 551
 552    -----------------
 553    -- Eigenvalues --
 554    -----------------
 555
 556    procedure Eigenvalues
 557      (A      : Real_Matrix;
 558       Values : out Real_Vector)
 559    is
 560       N      : constant Natural := Length (A);
 561       L_Work : Real_Vector (1 .. 1);
 562       Info   : aliased Integer;
 563
 564       B   : Real_Matrix (1 .. N, 1 .. N);
 565       Tau : Real_Vector (1 .. N);
 566       E   : Real_Vector (1 .. N);
 567       pragma Warnings (Off, B);
 568       pragma Warnings (Off, Tau);
 569       pragma Warnings (Off, E);
 570
 571    begin
 572       if Values'Length /= N then
 573          raise Constraint_Error with "wrong length for output vector";
 574       end if;
 575
 576       if N = 0 then
 577          return;
 578       end if;
 579
 580       --  Initialize working matrix and check for symmetric input matrix
 581
 582       Transpose (A, B);
 583
 584       if A /= B then
 585          raise Argument_Error with "matrix not symmetric";
 586       end if;
 587
 588       --  Find size of work area
 589
 590       sytrd (Uplo   => Lower'Access,
 591              N      => N,
 592              A      => B,
 593              Ld_A   => N,
 594              D      => Values,
 595              E      => E,
 596              Tau    => Tau,
 597              Work   => L_Work,
 598              L_Work => -1,
 599              Info   => Info'Access);
 600
 601       declare
 602          Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
 603          pragma Warnings (Off, Work);
 604
 605       begin
 606          --  Reduce matrix to tridiagonal form
 607
 608          sytrd (Uplo   => Lower'Access,
 609                 N      => N,
 610                 A      => B,
 611                 Ld_A   => A'Length (1),
 612                 D      => Values,
 613                 E      => E,
 614                 Tau    => Tau,
 615                 Work   => Work,
 616                 L_Work => Work'Length,
 617                 Info   => Info'Access);
 618
 619          if Info /= 0 then
 620             raise Constraint_Error;
 621          end if;
 622
 623          --  Compute all eigenvalues using QR algorithm
 624
 625          sterf (N      => N,
 626                 D      => Values,
 627                 E      => E,
 628                 Info   => Info'Access);
 629
 630          if Info /= 0 then
 631             raise Constraint_Error with
 632                "eigenvalues computation failed to converge";
 633          end if;
 634       end;
 635    end Eigenvalues;
 636
 637    function Eigenvalues (A : Real_Matrix) return Real_Vector is
 638       R : Real_Vector (A'Range (1));
 639    begin
 640       Eigenvalues (A, R);
 641       return R;
 642    end Eigenvalues;
 643
 644    -------------
 645    -- Inverse --
 646    -------------
 647
 648    procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is
 649       N      : constant Integer := Length (A);
 650       Piv    : Integer_Vector (1 .. N);
 651       L_Work : Real_Vector (1 .. 1);
 652       Info   : aliased Integer := -1;
 653
 654    begin
 655       --  All computations are done using column-major order, but this works
 656       --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
 657
 658       R := A;
 659
 660       --  Compute LU decomposition
 661
 662       getrf (M      => N,
 663              N      => N,
 664              A      => R,
 665              Ld_A   => N,
 666              I_Piv  => Piv,
 667              Info   => Info'Access);
 668
 669       if Info /= 0 then
 670          raise Constraint_Error with "inverting singular matrix";
 671       end if;
 672
 673       --  Determine size of work area
 674
 675       getri (N      => N,
 676              A      => R,
 677              Ld_A   => N,
 678              I_Piv  => Piv,
 679              Work   => L_Work,
 680              L_Work => -1,
 681              Info   => Info'Access);
 682
 683       if Info /= 0 then
 684          raise Constraint_Error;
 685       end if;
 686
 687       declare
 688          Work : Real_Vector (1 .. Integer (L_Work (1)));
 689          pragma Warnings (Off, Work);
 690
 691       begin
 692          --  Compute inverse from LU decomposition
 693
 694          getri (N      => N,
 695                 A      => R,
 696                 Ld_A   => N,
 697                 I_Piv  => Piv,
 698                 Work   => Work,
 699                 L_Work => Work'Length,
 700                 Info   => Info'Access);
 701
 702          if Info /= 0 then
 703             raise Constraint_Error with "inverting singular matrix";
 704          end if;
 705
 706          --  ??? Should iterate with gerfs, based on implementation advice
 707       end;
 708    end Inverse;
 709
 710    function Inverse (A : Real_Matrix) return Real_Matrix is
 711       R : Real_Matrix (A'Range (2), A'Range (1));
 712    begin
 713       Inverse (A, R);
 714       return R;
 715    end Inverse;
 716
 717    -----------
 718    -- Solve --
 719    -----------
 720
 721    procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector) is
 722    begin
 723       if Length (A) /= X'Length then
 724          raise Constraint_Error with
 725            "incompatible matrix and vector 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    procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix) is
 734    begin
 735       if Length (A) /= X'Length (1) then
 736          raise Constraint_Error with "incompatible matrix dimensions";
 737       end if;
 738
 739       --  ??? Should solve directly, is faster and more accurate
 740
 741       B := Inverse (A) * X;
 742    end Solve;
 743
 744    function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
 745       B : Real_Vector (A'Range (2));
 746    begin
 747       Solve (A, X, B);
 748       return B;
 749    end Solve;
 750
 751    function Solve (A, X : Real_Matrix) return Real_Matrix is
 752       B : Real_Matrix (A'Range (2), X'Range (2));
 753    begin
 754       Solve (A, X, B);
 755       return B;
 756    end Solve;
 757
 758    ---------------
 759    -- Transpose --
 760    ---------------
 761
 762    function Transpose (X : Real_Matrix) return Real_Matrix is
 763       R : Real_Matrix (X'Range (2), X'Range (1));
 764    begin
 765       Transpose (X, R);
 766
 767       return R;
 768    end Transpose;
 769
 770    -----------------
 771    -- Unit_Matrix --
 772    -----------------
 773
 774    function Unit_Matrix
 775      (Order   : Positive;
 776       First_1 : Integer := 1;
 777       First_2 : Integer := 1) return Real_Matrix
 778      renames Instantiations.Unit_Matrix;
 779
 780    -----------------
 781    -- Unit_Vector --
 782    -----------------
 783
 784    function Unit_Vector
 785      (Index : Integer;
 786       Order : Positive;
 787       First : Integer := 1) return Real_Vector
 788      renames Instantiations.Unit_Vector;
 789
 790 end Ada.Numerics.Generic_Real_Arrays;