2011-10-13 Geert Bosch <bosch@adacore.com>

	* a-ngrear.adb (Solve): Make generic and move to
	System.Generic_Array_Operations.
	* s-gearop.ads (Matrix_Vector_Solution, Matrix_Matrix_Solution):
	New generic solvers to	compute a vector resp. matrix Y such
	that A * Y = X, approximately.
	* s-gearop.adb (Matrix_Vector_Solution, Matrix_Matrix_Solution):
	Implement using Forward_Eliminate and Back_Substitute
	* a-ngcoar.adb: Reimplement in pure Ada to remove dependencies
	on BLAS and LAPACK.
	* a-ngcoar.ads ("abs"): Fix return type to be real.



git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@179912 138bc75d-0d04-0410-961f-82ee72b054a4

6 files changed, 293 insertions, 472 deletions

diff --git a/gcc/ada/ChangeLog b/gcc/ada/ChangeLog
index 461e4d1eaae..04b04fd7e8f 100644
--- a/gcc/ada/ChangeLog
+++ b/gcc/ada/ChangeLog
@@ -1,3 +1,16 @@
+2011-10-13  Geert Bosch  <bosch@adacore.com>
+
+	* a-ngrear.adb (Solve): Make generic and move to
+	System.Generic_Array_Operations.
+	* s-gearop.ads (Matrix_Vector_Solution, Matrix_Matrix_Solution):
+	New generic solvers to	compute a vector resp. matrix Y such
+	that A * Y = X, approximately.
+	* s-gearop.adb (Matrix_Vector_Solution, Matrix_Matrix_Solution):
+	Implement using Forward_Eliminate and Back_Substitute
+	* a-ngcoar.adb: Reimplement in pure Ada to remove dependencies
+	on BLAS and LAPACK.
+	* a-ngcoar.ads ("abs"): Fix return type to be real.
+
 2011-10-13  Eric Botcazou  <ebotcazou@adacore.com>
 
 	PR ada/50589
diff --git a/gcc/ada/a-ngcoar.adb b/gcc/ada/a-ngcoar.adb
index 9979d6baec6..3be88491d7f 100644
--- a/gcc/ada/a-ngcoar.adb
+++ b/gcc/ada/a-ngcoar.adb
@@ -6,7 +6,7 @@
 --                                                                          --
 --                                 B o d y                                  --
 --                                                                          --
---            Copyright (C) 2006-2009, Free Software Foundation, Inc.       --
+--            Copyright (C) 2006-2011, Free Software Foundation, Inc.       --
 --                                                                          --
 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
 -- terms of the  GNU General Public License as published  by the Free Soft- --
@@ -30,66 +30,35 @@
 ------------------------------------------------------------------------------
 
 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
-with System.Generic_Complex_BLAS;
-with System.Generic_Complex_LAPACK;
+with Ada.Numerics; use Ada.Numerics;
 
 package body Ada.Numerics.Generic_Complex_Arrays is
 
-   --  Operations involving inner products use BLAS library implementations.
-   --  This allows larger matrices and vectors to be computed efficiently,
-   --  taking into account memory hierarchy issues and vector instructions
-   --  that vary widely between machines.
-
    --  Operations that are defined in terms of operations on the type Real,
    --  such as addition, subtraction and scaling, are computed in the canonical
    --  way looping over all elements.
 
-   --  Operations for solving linear systems and computing determinant,
-   --  eigenvalues, eigensystem and inverse, are implemented using the
-   --  LAPACK library.
-
-   type BLAS_Real_Vector is array (Integer range <>) of Real;
-
-   package BLAS is new System.Generic_Complex_BLAS
-     (Real           => Real,
-      Complex_Types  => Complex_Types,
-      Complex_Vector => Complex_Vector,
-      Complex_Matrix => Complex_Matrix);
-
-   package LAPACK is new System.Generic_Complex_LAPACK
-     (Real           => Real,
-      Real_Vector    => BLAS_Real_Vector,
-      Complex_Types  => Complex_Types,
-      Complex_Vector => Complex_Vector,
-      Complex_Matrix => Complex_Matrix);
+   package Ops renames System.Generic_Array_Operations;
 
    subtype Real is Real_Arrays.Real;
    --  Work around visibility bug ???
 
-   use BLAS, LAPACK;
-
-   --  Procedure versions of functions returning unconstrained values.
-   --  This allows for inlining the function wrapper.
+   function Is_Non_Zero (X : Complex) return Boolean is (X /= (0.0, 0.0));
+   --  Needed by Back_Substitute
 
-   procedure Eigenvalues
-     (A      : Complex_Matrix;
-      Values : out Real_Vector);
+   procedure Back_Substitute is new Ops.Back_Substitute
+     (Scalar        => Complex,
+      Matrix        => Complex_Matrix,
+      Is_Non_Zero   => Is_Non_Zero);
 
-   procedure Inverse
-     (A      : Complex_Matrix;
-      R      : out Complex_Matrix);
+   procedure Forward_Eliminate is new Ops.Forward_Eliminate
+    (Scalar        => Complex,
+     Real          => Real'Base,
+     Matrix        => Complex_Matrix,
+     Zero          => (0.0, 0.0),
+     One           => (1.0, 0.0));
 
-   procedure Solve
-     (A      : Complex_Matrix;
-      X      : Complex_Vector;
-      B      : out Complex_Vector);
-
-   procedure Solve
-     (A      : Complex_Matrix;
-      X      : Complex_Matrix;
-      B      : out Complex_Matrix);
-
-   procedure Transpose is new System.Generic_Array_Operations.Transpose
+   procedure Transpose is new Ops.Transpose
                                 (Scalar => Complex,
                                  Matrix => Complex_Matrix);
 
@@ -98,6 +67,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is
 
    function Length is new Square_Matrix_Length (Complex, Complex_Matrix);
 
+   --  Instant a generic square root implementation here, in order to avoid
+   --  instantiating a complete copy of Generic_Elementary_Functions.
+   --  Speed of the square root is not a big concern here.
+
+   function Sqrt is new Ops.Sqrt (Real'Base);
+
    --  Instantiating the following subprograms directly would lead to
    --  name clashes, so use a local package.
 
@@ -155,6 +130,14 @@ package body Ada.Numerics.Generic_Complex_Arrays is
                              Right_Vector  => Complex_Vector,
                              Zero          => (0.0, 0.0));
 
+      function "*" is new Inner_Product
+                            (Left_Scalar   => Complex,
+                             Right_Scalar  => Complex,
+                             Result_Scalar => Complex,
+                             Left_Vector   => Complex_Vector,
+                             Right_Vector  => Complex_Vector,
+                             Zero          => (0.0, 0.0));
+
       function "*" is new Outer_Product
                             (Left_Scalar   => Complex,
                              Right_Scalar  => Complex,
@@ -229,6 +212,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is
                              Result_Vector => Complex_Vector,
                              Zero          => (0.0, 0.0));
 
+      function "*" is new Matrix_Vector_Product
+                            (Left_Scalar   => Complex,
+                             Right_Scalar  => Complex,
+                             Result_Scalar => Complex,
+                             Matrix        => Complex_Matrix,
+                             Right_Vector  => Complex_Vector,
+                             Result_Vector => Complex_Vector,
+                             Zero          => (0.0, 0.0));
+
       function "*" is new Vector_Matrix_Product
                             (Left_Scalar   => Real'Base,
                              Right_Scalar  => Complex,
@@ -247,6 +239,24 @@ package body Ada.Numerics.Generic_Complex_Arrays is
                              Result_Vector => Complex_Vector,
                              Zero          => (0.0, 0.0));
 
+      function "*" is new Vector_Matrix_Product
+                            (Left_Scalar   => Complex,
+                             Right_Scalar  => Complex,
+                             Result_Scalar => Complex,
+                             Left_Vector   => Complex_Vector,
+                             Matrix        => Complex_Matrix,
+                             Result_Vector => Complex_Vector,
+                             Zero          => (0.0, 0.0));
+
+      function "*" is new Matrix_Matrix_Product
+                            (Left_Scalar   => Complex,
+                             Right_Scalar  => Complex,
+                             Result_Scalar => Complex,
+                             Left_Matrix   => Complex_Matrix,
+                             Right_Matrix  => Complex_Matrix,
+                             Result_Matrix => Complex_Matrix,
+                             Zero          => (0.0, 0.0));
+
       function "*" is new Matrix_Matrix_Product
                             (Left_Scalar   => Real'Base,
                              Right_Scalar  => Complex,
@@ -445,6 +455,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is
                              Result_Matrix => Complex_Matrix,
                              Operation     => "/");
 
+      -----------
+      -- "abs" --
+      -----------
+
+      function "abs" is new L2_Norm
+                              (X_Scalar      => Complex,
+                               Result_Real   => Real'Base,
+                               X_Vector      => Complex_Vector);
+
       --------------
       -- Argument --
       --------------
@@ -671,6 +690,16 @@ package body Ada.Numerics.Generic_Complex_Arrays is
                              Y_Matrix      => Real_Matrix,
                              Update        => Set_Re);
 
+      -----------
+      -- Solve --
+      -----------
+
+      function Solve is
+         new Matrix_Vector_Solution (Complex, Complex_Vector, Complex_Matrix);
+
+      function Solve is
+         new Matrix_Matrix_Solution (Complex, Complex_Matrix);
+
       -----------------
       -- Unit_Matrix --
       -----------------
@@ -686,7 +715,6 @@ package body Ada.Numerics.Generic_Complex_Arrays is
                              Vector        => Complex_Vector,
                              Zero          => (0.0, 0.0),
                              One           => (1.0, 0.0));
-
    end Instantiations;
 
    ---------
@@ -696,15 +724,7 @@ package body Ada.Numerics.Generic_Complex_Arrays is
    function "*"
      (Left  : Complex_Vector;
       Right : Complex_Vector) return Complex
-   is
-   begin
-      if Left'Length /= Right'Length then
-         raise Constraint_Error with
-            "vectors are of different length in inner product";
-      end if;
-
-      return dot (Left'Length, X => Left, Y => Right);
-   end "*";
+     renames Instantiations."*";
 
    function "*"
      (Left  : Real_Vector;
@@ -738,31 +758,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
 
    function "*"
      (Left  : Complex_Matrix;
-      Right : Complex_Matrix)
-      return  Complex_Matrix
-   is
-      R : Complex_Matrix (Left'Range (1), Right'Range (2));
-
-   begin
-      if Left'Length (2) /= Right'Length (1) then
-         raise Constraint_Error with
-            "incompatible dimensions in matrix-matrix multiplication";
-      end if;
-
-      gemm (Trans_A => No_Trans'Access,
-            Trans_B => No_Trans'Access,
-            M       => Right'Length (2),
-            N       => Left'Length (1),
-            K       => Right'Length (1),
-            A       => Right,
-            Ld_A    => Right'Length (2),
-            B       => Left,
-            Ld_B    => Left'Length (2),
-            C       => R,
-            Ld_C    => R'Length (2));
-
-      return R;
-   end "*";
+      Right : Complex_Matrix) return  Complex_Matrix
+     renames Instantiations."*";
 
    function "*"
      (Left  : Complex_Vector;
@@ -772,48 +769,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is
    function "*"
      (Left  : Complex_Vector;
       Right : Complex_Matrix) return Complex_Vector
-   is
-      R : Complex_Vector (Right'Range (2));
-
-   begin
-      if Left'Length /= Right'Length (1) then
-         raise Constraint_Error with
-           "incompatible dimensions in vector-matrix multiplication";
-      end if;
-
-      gemv (Trans => No_Trans'Access,
-            M     => Right'Length (2),
-            N     => Right'Length (1),
-            A     => Right,
-            Ld_A  => Right'Length (2),
-            X     => Left,
-            Y     => R);
-
-      return R;
-   end "*";
+     renames Instantiations."*";
 
    function "*"
      (Left  : Complex_Matrix;
       Right : Complex_Vector) return Complex_Vector
-   is
-      R : Complex_Vector (Left'Range (1));
-
-   begin
-      if Left'Length (2) /= Right'Length then
-         raise Constraint_Error with
-            "incompatible dimensions in matrix-vector multiplication";
-      end if;
-
-      gemv (Trans => Trans'Access,
-            M     => Left'Length (2),
-            N     => Left'Length (1),
-            A     => Left,
-            Ld_A  => Left'Length (2),
-            X     => Right,
-            Y     => R);
-
-      return R;
-   end "*";
+     renames Instantiations."*";
 
    function "*"
      (Left  : Real_Matrix;
@@ -984,10 +945,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
    -- "abs" --
    -----------
 
-   function "abs" (Right : Complex_Vector) return Complex is
-   begin
-      return (nrm2 (Right'Length, Right), 0.0);
-   end "abs";
+   function "abs" (Right : Complex_Vector) return Real'Base
+      renames Instantiations."abs";
 
    --------------
    -- Argument --
@@ -1070,38 +1029,12 @@ package body Ada.Numerics.Generic_Complex_Arrays is
    -----------------
 
    function Determinant (A : Complex_Matrix) return Complex is
-      N    : constant Integer := Length (A);
-      LU   : Complex_Matrix (1 .. N, 1 .. N) := A;
-      Piv  : Integer_Vector (1 .. N);
-      Info : aliased Integer := -1;
-      Neg  : Boolean;
-      Det  : Complex;
-
+      M : Complex_Matrix := A;
+      B : Complex_Matrix (A'Range (1), 1 .. 0);
+      R : Complex;
    begin
-      if N = 0 then
-         return (0.0, 0.0);
-      end if;
-
-      getrf (N, N, LU, N, Piv, Info'Access);
-
-      if Info /= 0 then
-         raise Constraint_Error with "ill-conditioned matrix";
-      end if;
-
-      Det := LU (1, 1);
-      Neg := Piv (1) /= 1;
-
-      for J in 2 .. N loop
-         Det := Det * LU (J, J);
-         Neg := Neg xor (Piv (J) /= J);
-      end loop;
-
-      if Neg then
-         return -Det;
-
-      else
-         return Det;
-      end if;
+      Forward_Eliminate (M, B, R);
+      return R;
    end Determinant;
 
    -----------------
@@ -1113,174 +1046,96 @@ package body Ada.Numerics.Generic_Complex_Arrays is
       Values  : out Real_Vector;
       Vectors : out Complex_Matrix)
    is
-      Job_Z    : aliased Character := 'V';
-      Rng      : aliased Character := 'A';
-      Uplo     : aliased Character := 'U';
-
-      N        : constant Natural := Length (A);
-      W        : BLAS_Real_Vector (Values'Range);
-      M        : Integer;
-      B        : Complex_Matrix (1 .. N, 1 .. N);
-      L_Work   : Complex_Vector (1 .. 1);
-      LR_Work  : BLAS_Real_Vector (1 .. 1);
-      LI_Work  : Integer_Vector (1 .. 1);
-      I_Supp_Z : Integer_Vector (1 .. 2 * N);
-      Info     : aliased Integer;
+      N : constant Natural := Length (A);
+
+      --  For a Hermitian matrix C, we convert the eigenvalue problem to a
+      --  real symmetric one: if C = A + i * B, then the (N, N) complex
+      --  eigenvalue problem:
+      --     (A + i * B) * (u + i * v) = Lambda * (u + i * v)
+      --
+      --  is equivalent to the (2 * N, 2 * N) real eigenvalue problem:
+      --     [  A, B ] [ u ] = Lambda * [ u ]
+      --     [ -B, A ] [ v ]            [ v ]
+      --
+      --  Note that the (2 * N, 2 * N) matrix above is symmetric, as
+      --  Transpose (A) = A and Transpose (B) = -B if C is Hermitian.
+
+      --  We solve this eigensystem using the real-valued algorithms. The final
+      --  result will have every eigenvalue twice, so in the sorted output we
+      --  just pick every second value, with associated eigenvector u + i * v.
+
+      M    : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
+      Vals : Real_Vector (1 .. 2 * N);
+      Vecs : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
 
    begin
-      if Values'Length /= N then
-         raise Constraint_Error with "wrong length for output vector";
-      end if;
-
-      if Vectors'First (1) /= A'First (1)
-        or else Vectors'Last (1) /= A'Last (1)
-        or else Vectors'First (2) /= A'First (2)
-        or else Vectors'Last (2) /= A'Last (2)
-      then
-         raise Constraint_Error with "wrong dimensions for output matrix";
-      end if;
-
-      if N = 0 then
-         return;
-      end if;
-
-      --  Check for hermitian matrix ???
-      --  Copy only required triangle ???
-
-      B := A;
-
-      --  Find size of work area
-
-      heevr
-        (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
-         M        => M,
-         W        => W,
-         Z        => Vectors,
-         Ld_Z     => N,
-         I_Supp_Z => I_Supp_Z,
-         Work     => L_Work,
-         L_Work   => -1,
-         R_Work   => LR_Work,
-         LR_Work  => -1,
-         I_Work   => LI_Work,
-         LI_Work  => -1,
-         Info     => Info'Access);
-
-      if Info /= 0 then
-         raise Constraint_Error;
-      end if;
-
-      declare
-         Work   : Complex_Vector (1 .. Integer (L_Work (1).Re));
-         R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
-         I_Work : Integer_Vector (1 .. LI_Work (1));
-
-      begin
-         heevr
-           (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
-            M        => M,
-            W        => W,
-            Z        => Vectors,
-            Ld_Z     => N,
-            I_Supp_Z => I_Supp_Z,
-            Work     => Work,
-            L_Work   => Work'Length,
-            R_Work   => R_Work,
-            LR_Work  => LR_Work'Length,
-            I_Work   => I_Work,
-            LI_Work  => LI_Work'Length,
-            Info     => Info'Access);
-
-         if Info /= 0 then
-            raise Constraint_Error with "inverting non-Hermitian matrix";
-         end if;
-
-         for J in Values'Range loop
-            Values (J) := W (J);
+      for J in 1 .. N loop
+         for K in 1 .. N loop
+            declare
+               C : constant Complex :=
+                     (A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
+            begin
+               M (J, K) := Re (C);
+               M (J + N, K + N) := Re (C);
+               M (J + N, K) := Im (C);
+               M (J, K + N) := -Im (C);
+            end;
          end loop;
-      end;
+      end loop;
+
+      Eigensystem (M, Vals, Vecs);
+
+      for J in 1 .. N loop
+         declare
+            Col : constant Integer := Values'First + (J - 1);
+         begin
+            Values (Col) := Vals (2 * J);
+
+            for K in 1 .. N loop
+               declare
+                  Row : constant Integer := Vectors'First (2) + (K - 1);
+               begin
+                  Vectors (Row, Col)
+                     := (Vecs (J * 2, Col), Vecs (J * 2, Col + N));
+               end;
+            end loop;
+         end;
+      end loop;
    end Eigensystem;
 
    -----------------
    -- Eigenvalues --
    -----------------
 
-   procedure Eigenvalues
-     (A      : Complex_Matrix;
-      Values : out Real_Vector)
-   is
-      Job_Z    : aliased Character := 'N';
-      Rng      : aliased Character := 'A';
-      Uplo     : aliased Character := 'U';
-      N        : constant Natural := Length (A);
-      B        : Complex_Matrix (1 .. N, 1 .. N) := A;
-      Z        : Complex_Matrix (1 .. 1, 1 .. 1);
-      W        : BLAS_Real_Vector (Values'Range);
-      L_Work   : Complex_Vector (1 .. 1);
-      LR_Work  : BLAS_Real_Vector (1 .. 1);
-      LI_Work  : Integer_Vector (1 .. 1);
-      I_Supp_Z : Integer_Vector (1 .. 2 * N);
-      M        : Integer;
-      Info     : aliased Integer;
+   function Eigenvalues (A : Complex_Matrix) return Real_Vector is
+      --  See Eigensystem for a description of the algorithm
 
+      N : constant Natural := Length (A);
+      R : Real_Vector (A'Range (1));
+
+      M    : Real_Matrix (1 .. 2 * N, 1 .. 2 * N);
+      Vals : Real_Vector (1 .. 2 * N);
    begin
-      if Values'Length /= N then
-         raise Constraint_Error with "wrong length for output vector";
-      end if;
-
-      if N = 0 then
-         return;
-      end if;
-
-      --  Check for hermitian matrix ???
-
-      --  Find size of work area
-
-      heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
-             M        => M,
-             W        => W,
-             Z        => Z,
-             Ld_Z     => 1,
-             I_Supp_Z => I_Supp_Z,
-             Work     => L_Work,  L_Work  => -1,
-             R_Work   => LR_Work, LR_Work => -1,
-             I_Work   => LI_Work, LI_Work => -1,
-             Info     => Info'Access);
-
-      if Info /= 0 then
-         raise Constraint_Error;
-      end if;
-
-      declare
-         Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
-         R_Work : BLAS_Real_Vector (1 .. Integer (LR_Work (1)));
-         I_Work : Integer_Vector (1 .. LI_Work (1));
-      begin
-         heevr (Job_Z'Access, Rng'Access, Uplo'Access, N, B, N,
-                M        => M,
-                W        => W,
-                Z        => Z,
-                Ld_Z     => 1,
-                I_Supp_Z => I_Supp_Z,
-                Work     => Work,   L_Work  => Work'Length,
-                R_Work   => R_Work, LR_Work => R_Work'Length,
-                I_Work   => I_Work, LI_Work => I_Work'Length,
-                Info     => Info'Access);
-
-         if Info /= 0 then
-            raise Constraint_Error with "inverting singular matrix";
-         end if;
-
-         for J in Values'Range loop
-            Values (J) := W (J);
+      for J in 1 .. N loop
+         for K in 1 .. N loop
+            declare
+               C : constant Complex :=
+                     (A (A'First (1) + (J - 1), A'First (2) + (K - 1)));
+            begin
+               M (J, K) := Re (C);
+               M (J + N, K + N) := Re (C);
+               M (J + N, K) := Im (C);
+               M (J, K + N) := -Im (C);
+            end;
          end loop;
-      end;
-   end Eigenvalues;
+      end loop;
+
+      Vals := Eigenvalues (M);
+
+      for J in 1 .. N loop
+         R (A'First (1) + (J - 1)) := Vals (2 * J);
+      end loop;
 
-   function Eigenvalues (A : Complex_Matrix) return Real_Vector is
-      R : Real_Vector (A'Range (1));
-   begin
-      Eigenvalues (A, R);
       return R;
    end Eigenvalues;
 
@@ -1298,73 +1153,8 @@ package body Ada.Numerics.Generic_Complex_Arrays is
    -- Inverse --
    -------------
 
-   procedure Inverse (A : Complex_Matrix; R : out Complex_Matrix) is
-      N      : constant Integer := Length (A);
-      Piv    : Integer_Vector (1 .. N);
-      L_Work : Complex_Vector (1 .. 1);
-      Info   : aliased Integer := -1;
-
-   begin
-      --  All computations are done using column-major order, but this works
-      --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
-
-      R := A;
-
-      --  Compute LU decomposition
-
-      getrf (M      => N,
-             N      => N,
-             A      => R,
-             Ld_A   => N,
-             I_Piv  => Piv,
-             Info   => Info'Access);
-
-      if Info /= 0 then
-         raise Constraint_Error with "inverting singular matrix";
-      end if;
-
-      --  Determine size of work area
-
-      getri (N      => N,
-             A      => R,
-             Ld_A   => N,
-             I_Piv  => Piv,
-             Work   => L_Work,
-             L_Work => -1,
-             Info   => Info'Access);
-
-      if Info /= 0 then
-         raise Constraint_Error;
-      end if;
-
-      declare
-         Work : Complex_Vector (1 .. Integer (L_Work (1).Re));
-
-      begin
-         --  Compute inverse from LU decomposition
-
-         getri (N      => N,
-                A      => R,
-                Ld_A   => N,
-                I_Piv  => Piv,
-                Work   => Work,
-                L_Work => Work'Length,
-                Info   => Info'Access);
-
-         if Info /= 0 then
-            raise Constraint_Error with "inverting singular matrix";
-         end if;
-
-         --  ??? Should iterate with gerfs, based on implementation advice
-      end;
-   end Inverse;
-
    function Inverse (A : Complex_Matrix) return Complex_Matrix is
-      R : Complex_Matrix (A'Range (2), A'Range (1));
-   begin
-      Inverse (A, R);
-      return R;
-   end Inverse;
+     (Solve (A, Unit_Matrix (Length (A))));
 
    -------------
    -- Modulus --
@@ -1418,53 +1208,15 @@ package body Ada.Numerics.Generic_Complex_Arrays is
    -- Solve --
    -----------
 
-   procedure Solve
-     (A : Complex_Matrix;
-      X : Complex_Vector;
-      B : out Complex_Vector)
-   is
-   begin
-      if Length (A) /= X'Length then
-         raise Constraint_Error with
-           "incompatible matrix and vector dimensions";
-      end if;
-
-      --  ??? Should solve directly, is faster and more accurate
-
-      B := Inverse (A) * X;
-   end Solve;
-
-   procedure Solve
-     (A : Complex_Matrix;
-      X : Complex_Matrix;
-      B : out Complex_Matrix)
-   is
-   begin
-      if Length (A) /= X'Length (1) then
-         raise Constraint_Error with "incompatible matrix dimensions";
-      end if;
-
-      --  ??? Should solve directly, is faster and more accurate
-
-      B := Inverse (A) * X;
-   end Solve;
-
    function Solve
      (A : Complex_Matrix;
       X : Complex_Vector) return Complex_Vector
-   is
-      B : Complex_Vector (A'Range (2));
-   begin
-      Solve (A, X, B);
-      return B;
-   end Solve;
+     renames Instantiations.Solve;
 
-   function Solve (A, X : Complex_Matrix) return Complex_Matrix is
-      B : Complex_Matrix (A'Range (2), X'Range (2));
-   begin
-      Solve (A, X, B);
-      return B;
-   end Solve;
+   function Solve
+     (A : Complex_Matrix;
+      X : Complex_Matrix) return Complex_Matrix
+     renames Instantiations.Solve;
 
    ---------------
    -- Transpose --
diff --git a/gcc/ada/a-ngcoar.ads b/gcc/ada/a-ngcoar.ads
index abffbd1b636..8f8f37a7906 100644
--- a/gcc/ada/a-ngcoar.ads
+++ b/gcc/ada/a-ngcoar.ads
@@ -66,7 +66,7 @@ package Ada.Numerics.Generic_Complex_Arrays is
    function "+" (Left, Right : Complex_Vector) return Complex_Vector;
    function "-" (Left, Right : Complex_Vector) return Complex_Vector;
    function "*" (Left, Right : Complex_Vector) return Complex;
-   function "abs" (Right : Complex_Vector) return Complex;
+   function "abs" (Right : Complex_Vector) return Real'Base;
 
    --  Mixed Real_Vector and Complex_Vector arithmetic operations
 
diff --git a/gcc/ada/a-ngrear.adb b/gcc/ada/a-ngrear.adb
index c5ed66a3f7d..2a740b5c6b4 100644
--- a/gcc/ada/a-ngrear.adb
+++ b/gcc/ada/a-ngrear.adb
@@ -33,7 +33,7 @@
 --  reason for this is new Ada 2012 requirements that prohibit algorithms such
 --  as Strassen's algorithm, which may be used by some BLAS implementations. In
 --  addition, some platforms lacked suitable compilers to compile the reference
---  BLAS/LAPACK implementation. Finally, on some platforms there are be more
+--  BLAS/LAPACK implementation. Finally, on some platforms there are more
 --  floating point types than supported by BLAS/LAPACK.
 
 with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers;
@@ -337,6 +337,11 @@ package body Ada.Numerics.Generic_Real_Arrays is
            Result_Matrix => Real_Matrix,
            Operation     => "abs");
 
+      function Solve is
+         new Matrix_Vector_Solution (Real'Base, Real_Vector, Real_Matrix);
+
+      function Solve is new Matrix_Matrix_Solution (Real'Base, Real_Matrix);
+
       function Unit_Matrix is new
         Generic_Array_Operations.Unit_Matrix
           (Scalar        => Real'Base,
@@ -696,58 +701,11 @@ package body Ada.Numerics.Generic_Real_Arrays is
    -- Solve --
    -----------
 
-   function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
-      N   : constant Natural := Length (A);
-      MA  : Real_Matrix := A;
-      MX  : Real_Matrix (A'Range (1), 1 .. 1);
-      R   : Real_Vector (A'Range (2));
-      Det : Real'Base;
-
-   begin
-      if X'Length /= N then
-         raise Constraint_Error with "incompatible vector length";
-      end if;
-
-      for J in 0 .. MX'Length (1) - 1 loop
-         MX (MX'First (1) + J, 1) := X (X'First + J);
-      end loop;
-
-      Forward_Eliminate (MA, MX, Det);
-      Back_Substitute (MA, MX);
-
-      for J in 0 .. R'Length - 1 loop
-         R (R'First + J) := MX (MX'First (1) + J, 1);
-      end loop;
-
-      return R;
-   end Solve;
-
-   function Solve (A, X : Real_Matrix) return Real_Matrix is
-      N  : constant Natural := Length (A);
-      MA : Real_Matrix (A'Range (2), A'Range (2));
-      MB : Real_Matrix (A'Range (2), X'Range (2));
-      Det : Real'Base;
-
-   begin
-      if X'Length (1) /= N then
-         raise Constraint_Error with "matrices have unequal number of rows";
-      end if;
-
-      for J in 0 .. A'Length (1) - 1 loop
-         for K in MA'Range (2) loop
-            MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
-         end loop;
-
-         for K in MB'Range (2) loop
-            MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
-         end loop;
-      end loop;
-
-      Forward_Eliminate (MA, MB, Det);
-      Back_Substitute (MA, MB);
+   function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector
+      renames Instantiations.Solve;
 
-      return MB;
-   end Solve;
+   function Solve (A, X : Real_Matrix) return Real_Matrix
+      renames Instantiations.Solve;
 
    ----------------------
    -- Sort_Eigensystem --
diff --git a/gcc/ada/s-gearop.adb b/gcc/ada/s-gearop.adb
index 3aba5b9f450..58602e1e0a8 100644
--- a/gcc/ada/s-gearop.adb
+++ b/gcc/ada/s-gearop.adb
@@ -651,6 +651,75 @@ package body System.Generic_Array_Operations is
       return R;
    end  Matrix_Matrix_Product;
 
+   ----------------------------
+   -- Matrix_Vector_Solution --
+   ----------------------------
+
+   function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector is
+      N   : constant Natural := A'Length (1);
+      MA  : Matrix := A;
+      MX  : Matrix (A'Range (1), 1 .. 1);
+      R   : Vector (A'Range (2));
+      Det : Scalar;
+
+   begin
+      if A'Length (2) /= N then
+         raise Constraint_Error with "matrix is not square";
+      end if;
+
+      if X'Length /= N then
+         raise Constraint_Error with "incompatible vector length";
+      end if;
+
+      for J in 0 .. MX'Length (1) - 1 loop
+         MX (MX'First (1) + J, 1) := X (X'First + J);
+      end loop;
+
+      Forward_Eliminate (MA, MX, Det);
+      Back_Substitute (MA, MX);
+
+      for J in 0 .. R'Length - 1 loop
+         R (R'First + J) := MX (MX'First (1) + J, 1);
+      end loop;
+
+      return R;
+   end Matrix_Vector_Solution;
+
+   ----------------------------
+   -- Matrix_Matrix_Solution --
+   ----------------------------
+
+   function Matrix_Matrix_Solution (A, X : Matrix) return Matrix is
+      N  : constant Natural := A'Length (1);
+      MA : Matrix (A'Range (2), A'Range (2));
+      MB : Matrix (A'Range (2), X'Range (2));
+      Det : Scalar;
+
+   begin
+      if A'Length (2) /= N then
+         raise Constraint_Error with "matrix is not square";
+      end if;
+
+      if X'Length (1) /= N then
+         raise Constraint_Error with "matrices have unequal number of rows";
+      end if;
+
+      for J in 0 .. A'Length (1) - 1 loop
+         for K in MA'Range (2) loop
+            MA (MA'First (1) + J, K) := A (A'First (1) + J, K);
+         end loop;
+
+         for K in MB'Range (2) loop
+            MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
+         end loop;
+      end loop;
+
+      Forward_Eliminate (MA, MB, Det);
+      Back_Substitute (MA, MB);
+
+      return MB;
+   end Matrix_Matrix_Solution;
+
    ---------------------------
    -- Matrix_Vector_Product --
    ---------------------------
diff --git a/gcc/ada/s-gearop.ads b/gcc/ada/s-gearop.ads
index 9e9973c7d9c..f401da219e3 100644
--- a/gcc/ada/s-gearop.ads
+++ b/gcc/ada/s-gearop.ads
@@ -390,6 +390,35 @@ pragma Pure (Generic_Array_Operations);
      (Left  : Left_Matrix;
       Right : Right_Matrix) return Result_Matrix;
 
+   ----------------------------
+   -- Matrix_Vector_Solution --
+   ----------------------------
+
+   generic
+      type Scalar is private;
+      type Vector is array (Integer range <>) of Scalar;
+      type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+      with procedure Back_Substitute (M, N : in out Matrix) is <>;
+      with procedure Forward_Eliminate
+             (M   : in out Matrix;
+              N   : in out Matrix;
+              Det : out Scalar) is <>;
+   function Matrix_Vector_Solution (A : Matrix; X : Vector) return Vector;
+
+   ----------------------------
+   -- Matrix_Matrix_Solution --
+   ----------------------------
+
+   generic
+      type Scalar is private;
+      type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+      with procedure Back_Substitute (M, N : in out Matrix) is <>;
+      with procedure Forward_Eliminate
+             (M   : in out Matrix;
+              N   : in out Matrix;
+              Det : out Scalar) is <>;
+   function Matrix_Matrix_Solution (A : Matrix; X : Matrix) return Matrix;
+
    ----------
    -- Sqrt --
    ----------