Patchwork [Ada] Reimplement Ada.Numerics.Generic_Real_Arrays in pure Ada

login
register
mail settings
Submitter Arnaud Charlet
Date Aug. 29, 2011, 1:08 p.m.
Message ID <20110829130811.GA1834@adacore.com>
Download mbox | patch
Permalink /patch/112045/
State New
Headers show

Comments

Arnaud Charlet - Aug. 29, 2011, 1:08 p.m.
This new implementation avoids dependencies on BLAS and LAPACK.
The advantages are simplicity (no external requirements), more certainty
about accuracy and availability for all floating point types. Obvious
disadvantages are lack of target-specific optimizations that may affect
both accuracy, less flexible data formats and higher stack space
requirements and especially speed.

In particular, for solving eigensystems, Jacobi's method is used.
For larger matrices this may be one or two orders of magnitude
slower than more advanced algorithms, however, since it is expected
that these routines are generally used on smaller problems and because
this algorithm has the best stability and accuracy, this is an acceptable
trade-off.

Tested on x86_64-pc-linux-gnu, committed on trunk

2011-08-29  Geert Bosch  <bosch@adacore.com>

	* s-gearop.ads (Back_Substitute, Diagonal, Forward_Eliminate,
	L2_Norm, Swap_Column): New generic subprograms
	* s-gearop.adb (Back_Substitute, Diagonal, Forward_Eliminate,
	L2_Norm, Swap_Column): Implement new subprograms in order to
	eliminate dependency on BLAS and LAPACK libraries in
	Ada.Numerics.Generic_Real_Arrays and eventually also the complex
	version. Forward_Eliminate/Back_Substitute can be used to put a
	matrix in row echelon or reduced row echelon form using partial
	pivoting.
	* a-ngrear.adb: (Back_Substitute, Diagonal, Forward_Eleminate,
	Swap_Column): Instantiate from System.Generic_Array_Operations.
	("*", "abs"): Implement by instantiation from Generic_Array_Operations.
	(Sqrt): Local function for simple computation of square root without
	adding dependencies on Generic_Elementary_Functions.
	(Swap): New subprogram to exchange floating point numbers.
	(Inverse): Reimplement using Jordan-Gauss elimination.
	(Jacobi): New procedure implementing Jacobi's method for computation
	of eigensystems, based on Rutishauser's implementation.
	(L2_Norm): Implement directly using the inner product.
	(Sort_Eigensystem): Sort eigenvalue/eigenvector pairs in order of
	decreasing eigenvalue as required by the Ada RM.
	(Swap_Column): New helper procedure for Sort_Eigensystem.
	Remove with of System.Generic_Real_BLAS and System.Generic_Real_LAPACK.
	Add with of Ada.Containers.Generic_Anonymous_Array_Sort, for
	Sort_Eigensystems.
Pistola - Aug. 30, 2011, 8:25 a.m.
Dear Arnaud,


Arnaud Charlet wrote:
> 
> +--  This version of Generic_Real_Arrays avoids the use of BLAS and
> LAPACK. One
> +--  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
> 

May I ask why "Ada requirements prohibit algorithms such as Strassen's"?

Thanks,
mqk

Patch

Index: a-ngrear.adb
===================================================================
--- a-ngrear.adb	(revision 178155)
+++ a-ngrear.adb	(working copy)
@@ -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- --
@@ -29,52 +29,155 @@ 
 --                                                                          --
 ------------------------------------------------------------------------------
 
+--  This version of Generic_Real_Arrays avoids the use of BLAS and LAPACK. One
+--  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 many platforms there may be more
+--  floating point types than supported by BLAS/LAPACK.
+
+with Ada.Containers.Generic_Anonymous_Array_Sort; use Ada.Containers;
+
 with System; use System;
-with System.Generic_Real_BLAS;
-with System.Generic_Real_LAPACK;
 with System.Generic_Array_Operations; use System.Generic_Array_Operations;
 
 package body Ada.Numerics.Generic_Real_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.
+   package Ops renames System.Generic_Array_Operations;
 
-   --  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.
+   function Is_Non_Zero (X : Real'Base) return Boolean is (X /= 0.0);
 
-   --  Operations for solving linear systems and computing determinant,
-   --  eigenvalues, eigensystem and inverse, are implemented using the
-   --  LAPACK library.
+   procedure Back_Substitute is new Ops.Back_Substitute
+       (Scalar        => Real'Base,
+        Matrix        => Real_Matrix,
+        Is_Non_Zero   => Is_Non_Zero);
 
-   package BLAS is
-      new Generic_Real_BLAS (Real'Base, Real_Vector, Real_Matrix);
+   function Diagonal is new Ops.Diagonal
+        (Scalar       => Real'Base,
+         Vector       => Real_Vector,
+         Matrix       => Real_Matrix);
 
-   package LAPACK is
-      new Generic_Real_LAPACK (Real'Base, Real_Vector, Real_Matrix);
+   procedure Forward_Eliminate is new Ops.Forward_Eliminate
+       (Scalar        => Real'Base,
+        Matrix        => Real_Matrix,
+        Zero          => 0.0,
+        One           => 1.0);
 
-   use BLAS, LAPACK;
+   procedure Swap_Column is new Ops.Swap_Column
+       (Scalar        => Real'Base,
+        Matrix        => Real_Matrix);
 
-   --  Procedure versions of functions returning unconstrained values.
-   --  This allows for inlining the function wrapper.
-
-   procedure Eigenvalues (A : Real_Matrix; Values : out Real_Vector);
-   procedure Inverse   (A : Real_Matrix; R : out Real_Matrix);
-   procedure Solve     (A : Real_Matrix; X : Real_Vector; B : out Real_Vector);
-   procedure Solve     (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix);
-
-   procedure Transpose is new
-     Generic_Array_Operations.Transpose
+   procedure Transpose is new  Ops.Transpose
        (Scalar        => Real'Base,
         Matrix        => Real_Matrix);
 
+   function Is_Symmetric (A : Real_Matrix) return Boolean is
+     (Transpose (A) = A);
+   --  Return True iff A is symmetric, see RM G.3.1 (90).
+
+   function Is_Tiny (Value, Compared_To : Real) return Boolean is
+     (abs Compared_To + 100.0 * abs (Value) = abs Compared_To);
+   --  Return True iff the Value is much smaller in magnitude than the least
+   --  significant digit of Compared_To.
+
+   procedure Jacobi
+     (A               : Real_Matrix;
+      Values          : out Real_Vector;
+      Vectors         : out Real_Matrix;
+      Compute_Vectors : Boolean := True);
+   --  Perform Jacobi's eigensystem algorithm on real symmetric matrix A
+
+   function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
    --  Helper function that raises a Constraint_Error is the argument is
    --  not a square matrix, and otherwise returns its length.
 
-   function Length is new Square_Matrix_Length (Real'Base, Real_Matrix);
+   procedure Rotate (X, Y : in out Real; Sin, Tau : Real);
+   --  Perform a Givens rotation
 
+   procedure Sort_Eigensystem
+     (Values  : in out Real_Vector;
+      Vectors : in out Real_Matrix);
+   --  Sort Values and associated Vectors by decreasing absolute value
+
+   procedure Swap (Left, Right : in out Real);
+   --  Exchange Left and Right.
+
+   function Sqrt (X : Real) return Real;
+   --  Sqrt is implemented locally here, in order to avoid dragging in all of
+   --  the elementary functions. Speed of the square root is not a big concern
+   --  here. This also avoids depending on a specific floating point type.
+
+   ------------
+   -- Rotate --
+   ------------
+
+   procedure Rotate (X, Y : in out Real; Sin, Tau : Real) is
+      Old_X : constant Real := X;
+      Old_Y : constant Real := Y;
+   begin
+      X := Old_X - Sin * (Old_Y + Old_X * Tau);
+      Y := Old_Y + Sin * (Old_X - Old_Y * Tau);
+   end Rotate;
+
+   ----------
+   -- Sqrt --
+   ----------
+
+   function Sqrt (X : Real) return Real is
+      Root, Next : Real;
+
+   begin
+      --  Be defensive: any comparisons with NaN values will yield False.
+
+      if not (X > 0.0) then
+         if X = 0.0 then
+            return X;
+
+         else
+            raise Argument_Error;
+         end if;
+      end if;
+
+      --  Compute an initial estimate based on:
+
+      --     X = M * R**E and Sqrt (X) = Sqrt (M) * R**(E / 2.0),
+
+      --  where M is the mantissa, R is the radix and E the exponent.
+
+      --  By ignoring the mantissa and ignoring the case of an odd
+      --  exponent, we get a final error that is at most R. In other words,
+      --  the result has about a single bit precision.
+
+      Root := Real (Real'Machine_Radix) ** (Real'Exponent (X) / 2);
+
+      --  Because of the poor initial estimate, use the Babylonian method of
+      --  computing the square root, as it is stable for all inputs. Every step
+      --  will roughly double the precision of the result. Just a few steps
+      --  suffice in most cases. Eight iterations should give about 2**8 bits
+      --  of precision.
+
+      for J in 1 .. 8 loop
+         Next := (Root + X / Root) / 2.0;
+
+         exit when Root = Next;
+
+         Root := Next;
+      end loop;
+
+      return Root;
+   end Sqrt;
+
+   ----------
+   -- Swap --
+   ----------
+
+   procedure Swap (Left, Right : in out Real) is
+      Temp : constant Real := Left;
+   begin
+      Left := Right;
+      Right := Temp;
+   end Swap;
+
    --  Instantiating the following subprograms directly would lead to
    --  name clashes, so use a local package.
 
@@ -197,6 +300,45 @@ 
            Right_Vector  => Real_Vector,
            Matrix        => Real_Matrix);
 
+      function "*" is new
+        Inner_Product
+          (Left_Scalar   => Real'Base,
+           Right_Scalar  => Real'Base,
+           Result_Scalar => Real'Base,
+           Left_Vector   => Real_Vector,
+           Right_Vector  => Real_Vector,
+           Zero          => 0.0);
+
+      function "*" is new
+        Matrix_Vector_Product
+          (Left_Scalar   => Real'Base,
+           Right_Scalar  => Real'Base,
+           Result_Scalar => Real'Base,
+           Matrix        => Real_Matrix,
+           Right_Vector  => Real_Vector,
+           Result_Vector => Real_Vector,
+           Zero          => 0.0);
+
+      function "*" is new
+        Vector_Matrix_Product
+          (Left_Scalar   => Real'Base,
+           Right_Scalar  => Real'Base,
+           Result_Scalar => Real'Base,
+           Left_Vector   => Real_Vector,
+           Matrix        => Real_Matrix,
+           Result_Vector => Real_Vector,
+           Zero          => 0.0);
+
+      function "*" is new
+        Matrix_Matrix_Product
+          (Left_Scalar   => Real'Base,
+           Right_Scalar  => Real'Base,
+           Result_Scalar => Real'Base,
+           Left_Matrix   => Real_Matrix,
+           Right_Matrix  => Real_Matrix,
+           Result_Matrix => Real_Matrix,
+           Zero          => 0.0);
+
       function "/" is new
         Vector_Scalar_Elementwise_Operation
           (Left_Scalar   => Real'Base,
@@ -216,6 +358,13 @@ 
            Operation     => "/");
 
       function "abs" is new
+        L2_Norm
+          (Scalar        => Real'Base,
+           Vector        => Real_Vector,
+           Inner_Product => "*",
+           Sqrt          => Sqrt);
+
+      function "abs" is new
         Vector_Elementwise_Operation
           (X_Scalar      => Real'Base,
            Result_Scalar => Real'Base,
@@ -299,91 +448,23 @@ 
 
    --  Vector multiplication
 
-   function "*" (Left, Right : Real_Vector) return Real'Base is
-   begin
-      if Left'Length /= Right'Length then
-         raise Constraint_Error with
-            "vectors are of different length in inner product";
-      end if;
+   function "*" (Left, Right : Real_Vector) return Real'Base
+      renames Instantiations."*";
 
-      return dot (Left'Length, X => Left, Y => Right);
-   end "*";
-
    function "*" (Left, Right : Real_Vector) return Real_Matrix
       renames Instantiations."*";
 
-   function "*"
-     (Left : Real_Vector;
-      Right : Real_Matrix) return Real_Vector
-   is
-      R : Real_Vector (Right'Range (2));
+   function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector
+      renames Instantiations."*";
 
-   begin
-      if Left'Length /= Right'Length (1) then
-         raise Constraint_Error with
-           "incompatible dimensions in vector-matrix multiplication";
-      end if;
+   function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector
+      renames Instantiations."*";
 
-      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 "*";
-
-   function "*"
-     (Left : Real_Matrix;
-      Right : Real_Vector) return Real_Vector
-   is
-      R : Real_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 "*";
-
    --  Matrix Multiplication
 
-   function "*" (Left, Right : Real_Matrix) return Real_Matrix is
-      R : Real_Matrix (Left'Range (1), Right'Range (2));
+   function "*" (Left, Right : Real_Matrix) return Real_Matrix
+      renames Instantiations."*";
 
-   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 "*";
-
    ---------
    -- "/" --
    ---------
@@ -398,10 +479,8 @@ 
    -- "abs" --
    -----------
 
-   function "abs" (Right : Real_Vector) return Real'Base is
-   begin
-      return nrm2 (Right'Length, Right);
-   end "abs";
+   function "abs" (Right : Real_Vector) return Real'Base
+      renames Instantiations."abs";
 
    function "abs" (Right : Real_Vector) return Real_Vector
       renames Instantiations."abs";
@@ -414,29 +493,14 @@ 
    -----------------
 
    function Determinant (A : Real_Matrix) return Real'Base is
-      N    : constant Integer := Length (A);
-      LU   : Real_Matrix (1 .. N, 1 .. N) := A;
-      Piv  : Integer_Vector (1 .. N);
-      Info : aliased Integer := -1;
-      Det  : Real := 1.0;
+      M : Real_Matrix := A;
+      B : Real_Matrix (A'Range (1), 1 .. 0);
+      R : Real'Base;
 
    begin
-      getrf (M     => N,
-             N     => N,
-             A     => LU,
-             Ld_A  => N,
-             I_Piv => Piv,
-             Info  => Info'Access);
+      Forward_Eliminate (M, B, R);
 
-      if Info /= 0 then
-         raise Constraint_Error with "ill-conditioned matrix";
-      end if;
-
-      for J in 1 .. N loop
-         Det := (if Piv (J) /= J then -Det * LU (J, J) else Det * LU (J, J));
-      end loop;
-
-      return Det;
+      return R;
    end Determinant;
 
    -----------------
@@ -448,306 +512,319 @@ 
       Values  : out Real_Vector;
       Vectors : out Real_Matrix)
    is
-      N      : constant Natural := Length (A);
-      Tau    : Real_Vector (1 .. N);
-      L_Work : Real_Vector (1 .. 1);
-      Info   : aliased Integer;
+   begin
+      Jacobi (A, Values, Vectors, Compute_Vectors => True);
+      Sort_Eigensystem (Values, Vectors);
+   end Eigensystem;
 
-      E : Real_Vector (1 .. N);
-      pragma Warnings (Off, E);
+   -----------------
+   -- Eigenvalues --
+   -----------------
 
+   function Eigenvalues (A : Real_Matrix) return Real_Vector is
+      Values  : Real_Vector (A'Range (1));
+      Vectors : Real_Matrix (1 .. 0, 1 .. 0);
    begin
-      if Values'Length /= N then
-         raise Constraint_Error with "wrong length for output vector";
-      end if;
+      Jacobi (A, Values, Vectors, Compute_Vectors => False);
+      Sort_Eigensystem (Values, Vectors);
 
-      if N = 0 then
-         return;
-      end if;
+      return Values;
+   end Eigenvalues;
 
-      --  Initialize working matrix and check for symmetric input matrix
+   -------------
+   -- Inverse --
+   -------------
 
-      Transpose (A, Vectors);
+   function Inverse (A : Real_Matrix) return Real_Matrix is
+     (Solve (A, Unit_Matrix (Length (A))));
 
-      if A /= Vectors then
-         raise Argument_Error with "matrix not symmetric";
-      end if;
+   ------------
+   -- Jacobi --
+   ------------
 
-      --  Compute size of additional working space
+   procedure Jacobi
+     (A               : Real_Matrix;
+      Values          : out Real_Vector;
+      Vectors         : out Real_Matrix;
+      Compute_Vectors : Boolean := True)
+   is
+      --  This subprogram uses Carl Gustav Jacob Jacobi's iterative method
+      --  for computing eigenvalues and eigenvectors and is based on
+      --  Rutishauser's implementation.
 
-      sytrd (Uplo   => Lower'Access,
-             N      => N,
-             A      => Vectors,
-             Ld_A   => N,
-             D      => Values,
-             E      => E,
-             Tau    => Tau,
-             Work   => L_Work,
-             L_Work => -1,
-             Info   => Info'Access);
+      --  The given real symmetric matrix is transformed iteratively to
+      --  diagonal form through a sequence of appropriately chosen elementary
+      --  orthogonal transformations, called Jacobi rotations here.
 
-      declare
-         Work : Real_Vector (1 .. Integer'Max (Integer (L_Work (1)), 2 * N));
-         pragma Warnings (Off, Work);
+      --  The Jacobi method produces a systematic decrease of the sum of the
+      --  squares of off-diagonal elements. Convergence to zero is quadratic,
+      --  both for this implementation, as for the classic method that doesn't
+      --  use row-wise scanning for pivot selection.
 
-         Comp_Z : aliased constant Character := 'V';
+      --  The numerical stability and accuracy of Jacobi's method make it the
+      --  best choice here, even though for large matrices other methods will
+      --  be significantly more efficient in both time and space.
 
-      begin
-         --  Reduce matrix to tridiagonal form
+      --  While the eigensystem computations are absolutely foolproof for all
+      --  real symmetric matrices, in presence of invalid values, or similar
+      --  exceptional situations it might not. In such cases the results cannot
+      --  be trusted and Constraint_Error is raised.
 
-         sytrd (Uplo   => Lower'Access,
-                N      => N,
-                A      => Vectors,
-                Ld_A   => A'Length (1),
-                D      => Values,
-                E      => E,
-                Tau    => Tau,
-                Work   => Work,
-                L_Work => Work'Length,
-                Info   => Info'Access);
+      --  Note: this implementation needs temporary storage for 2 * N + N**2
+      --        values of type Real.
 
-         if Info /= 0 then
-            raise Program_Error;
-         end if;
+      Max_Iterations  : constant := 50;
 
-         --  Generate the real orthogonal matrix determined by sytrd
+      N               : constant Natural := Length (A);
 
-         orgtr (Uplo   => Lower'Access,
-                N      => N,
-                A      => Vectors,
-                Ld_A   => N,
-                Tau    => Tau,
-                Work   => Work,
-                L_Work => Work'Length,
-                Info   => Info'Access);
+      subtype Square_Matrix is Real_Matrix (1 .. N, 1 .. N);
 
-         if Info /= 0 then
-            raise Program_Error;
-         end if;
+      --  In order to annihilate the M (Row, Col) element, the
+      --  rotation parameters Cos and Sin are computed as
+      --  follows:
 
-         --  Compute all eigenvalues and eigenvectors using QR algorithm
+      --    Theta = Cot (2.0 * Phi)
+      --          = (Diag (Col) - Diag (Row)) / (2.0 * M (Row, Col))
 
-         steqr (Comp_Z => Comp_Z'Access,
-                N      => N,
-                D      => Values,
-                E      => E,
-                Z      => Vectors,
-                Ld_Z   => N,
-                Work   => Work,
-                Info   => Info'Access);
+      --  Then Tan (Phi) as the smaller root (in modulus) of
 
-         if Info /= 0 then
-            raise Constraint_Error with
-               "eigensystem computation failed to converge";
-         end if;
-      end;
-   end Eigensystem;
+      --    T**2 + 2 * T * Theta = 1 (or 0.5 / Theta, if Theta is large)
 
-   -----------------
-   -- Eigenvalues --
-   -----------------
+      function Compute_Tan (Theta : Real) return Real is
+         (Real'Copy_Sign (1.0 / (abs Theta + Sqrt (1.0 + Theta**2)), Theta));
 
-   procedure Eigenvalues
-     (A      : Real_Matrix;
-      Values : out Real_Vector)
-   is
-      N      : constant Natural := Length (A);
-      L_Work : Real_Vector (1 .. 1);
-      Info   : aliased Integer;
+      function Compute_Tan (P, H : Real) return Real is
+         (if Is_Tiny (P, Compared_To => H) then P / H
+          else Compute_Tan (Theta => H / (2.0 * P)));
 
-      B   : Real_Matrix (1 .. N, 1 .. N);
-      Tau : Real_Vector (1 .. N);
-      E   : Real_Vector (1 .. N);
-      pragma Warnings (Off, B);
-      pragma Warnings (Off, Tau);
-      pragma Warnings (Off, E);
+      function Sum_Strict_Upper (M : Square_Matrix) return Real;
+      --  Return the sum of all elements in the strict upper triangle of M
 
-   begin
-      if Values'Length /= N then
-         raise Constraint_Error with "wrong length for output vector";
-      end if;
+      ----------------------
+      -- Sum_Strict_Upper --
+      ----------------------
 
-      if N = 0 then
-         return;
-      end if;
+      function Sum_Strict_Upper (M : Square_Matrix) return Real is
+         Sum : Real := 0.0;
+      begin
+         for Row in 1 .. N - 1 loop
+            for Col in Row + 1 .. N loop
+               Sum := Sum + abs M (Row, Col);
+            end loop;
+         end loop;
 
-      --  Initialize working matrix and check for symmetric input matrix
+         return Sum;
+      end Sum_Strict_Upper;
 
-      Transpose (A, B);
+      M         : Square_Matrix := A; --  Work space for solving eigensystem
+      Threshold : Real;
+      Sum       : Real;
+      Diag      : Real_Vector (1 .. N);
+      Diag_Adj  : Real_Vector (1 .. N);
 
-      if A /= B then
-         raise Argument_Error with "matrix not symmetric";
-      end if;
+      --  The vector Diag_Adj indicates the amount of change in each value,
+      --  while Diag tracks the value itself and Values holds the values as
+      --  they were at the beginning. As the changes typically will be small
+      --  compared to the absolute value of Diag, at the end of each iteration
+      --  Diag is computed as Diag + Diag_Adj thus avoiding accumulating
+      --  rounding errors. This technique is due to Rutishauser.
 
-      --  Find size of work area
+   begin
+      if Compute_Vectors
+         and then (Vectors'Length (1) /= N or else Vectors'Length (2) /= N)
+      then
+         raise Constraint_Error with "incompatible matrix dimensions";
 
-      sytrd (Uplo   => Lower'Access,
-             N      => N,
-             A      => B,
-             Ld_A   => N,
-             D      => Values,
-             E      => E,
-             Tau    => Tau,
-             Work   => L_Work,
-             L_Work => -1,
-             Info   => Info'Access);
+      elsif Values'Length /= N then
+         raise Constraint_Error with "incompatible vector length";
 
-      declare
-         Work : Real_Vector (1 .. Integer'Min (Integer (L_Work (1)), 4 * N));
-         pragma Warnings (Off, Work);
+      elsif not Is_Symmetric (M) then
+         raise Constraint_Error with "matrix not symmetric";
+      end if;
 
-      begin
-         --  Reduce matrix to tridiagonal form
+      --  Note: Only the locally declared matrix M and vectors (Diag, Diag_Adj)
+      --        have lower bound equal to 1. The Vectors matrix may have
+      --        different bounds, so take care indexing elements. Assignment
+      --        as a whole is fine as sliding is automatic in that case.
 
-         sytrd (Uplo   => Lower'Access,
-                N      => N,
-                A      => B,
-                Ld_A   => A'Length (1),
-                D      => Values,
-                E      => E,
-                Tau    => Tau,
-                Work   => Work,
-                L_Work => Work'Length,
-                Info   => Info'Access);
+      Vectors := (if not Compute_Vectors then (1 .. 0 => (1 .. 0 => 0.0))
+                  else Unit_Matrix (Vectors'Length (1), Vectors'Length (2)));
+      Values := Diagonal (M);
 
-         if Info /= 0 then
-            raise Constraint_Error;
-         end if;
+      Sweep : for Iteration in 1 .. Max_Iterations loop
 
-         --  Compute all eigenvalues using QR algorithm
+         --  The first three iterations, perform rotation for any non-zero
+         --  element. After this, rotate only for those that are not much
+         --  smaller than the average off-diagnal element. After the fifth
+         --  iteration, additionally zero out off-diagonal elements that are
+         --  very small compared to elements on the diagonal with the same
+         --  column or row index.
 
-         sterf (N      => N,
-                D      => Values,
-                E      => E,
-                Info   => Info'Access);
+         Sum := Sum_Strict_Upper (M);
 
-         if Info /= 0 then
-            raise Constraint_Error with
-               "eigenvalues computation failed to converge";
-         end if;
-      end;
-   end Eigenvalues;
+         exit Sweep when Sum = 0.0;
 
-   function Eigenvalues (A : Real_Matrix) return Real_Vector is
-      R : Real_Vector (A'Range (1));
-   begin
-      Eigenvalues (A, R);
-      return R;
-   end Eigenvalues;
+         Threshold := (if Iteration < 4 then 0.2 * Sum / Real (N**2) else 0.0);
 
-   -------------
-   -- Inverse --
-   -------------
+         --  Iterate over all off-diagonal elements, rotating any that have
+         --  an absolute value that exceeds the threshold.
 
-   procedure Inverse (A : Real_Matrix; R : out Real_Matrix) is
-      N      : constant Integer := Length (A);
-      Piv    : Integer_Vector (1 .. N);
-      L_Work : Real_Vector (1 .. 1);
-      Info   : aliased Integer := -1;
+         Diag := Values;
+         Diag_Adj := (others => 0.0); -- Accumulates adjustments to Diag
 
-   begin
-      --  All computations are done using column-major order, but this works
-      --  out fine, because Transpose (Inverse (Transpose (A))) = Inverse (A).
+         for Row in 1 .. N - 1 loop
+            for Col in Row + 1 .. N loop
 
-      R := A;
+               --  If, before the rotation M (Row, Col) is tiny compared to
+               --  Diag (Row) and Diag (Col), rotation is skipped. This is
+               --  meaningful, as it produces no larger error than would be
+               --  produced anyhow if the rotation had been performed.
+               --  Suppress this optimization in the first four sweeps, so
+               --  that this procedure can be used for computing eigenvectors
+               --  of perturbed diagonal matrices.
 
-      --  Compute LU decomposition
+               if Iteration > 4
+                  and then Is_Tiny (M (Row, Col), Compared_To => Diag (Row))
+                  and then Is_Tiny (M (Row, Col), Compared_To => Diag (Col))
+               then
+                  M (Row, Col) := 0.0;
 
-      getrf (M      => N,
-             N      => N,
-             A      => R,
-             Ld_A   => N,
-             I_Piv  => Piv,
-             Info   => Info'Access);
+               elsif abs M (Row, Col) > Threshold then
+                  Perform_Rotation : declare
+                     Tan : constant Real := Compute_Tan (M (Row, Col),
+                                               Diag (Col) - Diag (Row));
+                     Cos : constant Real := 1.0 / Sqrt (1.0 + Tan**2);
+                     Sin : constant Real := Tan * Cos;
+                     Tau : constant Real := Sin / (1.0 + Cos);
+                     Adj : constant Real := Tan * M (Row, Col);
 
-      if Info /= 0 then
-         raise Constraint_Error with "inverting singular matrix";
-      end if;
+                  begin
+                     Diag_Adj (Row) := Diag_Adj (Row) - Adj;
+                     Diag_Adj (Col) := Diag_Adj (Col) + Adj;
+                     Diag (Row) := Diag (Row) - Adj;
+                     Diag (Col) := Diag (Col) + Adj;
 
-      --  Determine size of work area
+                     M (Row, Col) := 0.0;
 
-      getri (N      => N,
-             A      => R,
-             Ld_A   => N,
-             I_Piv  => Piv,
-             Work   => L_Work,
-             L_Work => -1,
-             Info   => Info'Access);
+                     for J in 1 .. Row - 1 loop        --  1 <= J < Row
+                        Rotate (M (J, Row), M (J, Col), Sin, Tau);
+                     end loop;
 
-      if Info /= 0 then
-         raise Constraint_Error;
-      end if;
+                     for J in Row + 1 .. Col - 1 loop  --  Row < J < Col
+                        Rotate (M (Row, J), M (J, Col), Sin, Tau);
+                     end loop;
 
-      declare
-         Work : Real_Vector (1 .. Integer (L_Work (1)));
-         pragma Warnings (Off, Work);
+                     for J in Col + 1 .. N loop        --  Col < J <= N
+                        Rotate (M (Row, J), M (Col, J), Sin, Tau);
+                     end loop;
 
-      begin
-         --  Compute inverse from LU decomposition
+                     for J in Vectors'Range (1) loop
+                        Rotate (Vectors (J, Row - 1 + Vectors'First (2)),
+                                Vectors (J, Col - 1 + Vectors'First (2)),
+                                Sin, Tau);
+                     end loop;
+                  end Perform_Rotation;
+               end if;
+            end loop;
+         end loop;
 
-         getri (N      => N,
-                A      => R,
-                Ld_A   => N,
-                I_Piv  => Piv,
-                Work   => Work,
-                L_Work => Work'Length,
-                Info   => Info'Access);
+         Values := Values + Diag_Adj;
+      end loop Sweep;
 
-         if Info /= 0 then
-            raise Constraint_Error with "inverting singular matrix";
-         end if;
+      --  All normal matrices with valid values should converge perfectly.
 
-         --  ??? Should iterate with gerfs, based on implementation advice
-      end;
-   end Inverse;
+      if Sum /= 0.0 then
+         raise Constraint_Error with "eigensystem solution does not converge";
+      end if;
+   end Jacobi;
 
-   function Inverse (A : Real_Matrix) return Real_Matrix is
-      R : Real_Matrix (A'Range (2), A'Range (1));
-   begin
-      Inverse (A, R);
-      return R;
-   end Inverse;
-
    -----------
    -- Solve --
    -----------
 
-   procedure Solve (A : Real_Matrix; X : Real_Vector; B : out Real_Vector) is
+   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 Length (A) /= X'Length then
-         raise Constraint_Error with
-           "incompatible matrix and vector dimensions";
+      if X'Length /= N then
+         raise Constraint_Error with "incompatible vector length";
       end if;
 
-      --  ??? Should solve directly, is faster and more accurate
+      for J in 0 .. MX'Length (1) - 1 loop
+         MX (MX'First (1) + J, 1) := X (X'First + J);
+      end loop;
 
-      B := Inverse (A) * X;
+      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;
 
-   procedure Solve (A : Real_Matrix; X : Real_Matrix; B : out Real_Matrix) is
+   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 Length (A) /= X'Length (1) then
-         raise Constraint_Error with "incompatible matrix dimensions";
+      if X'Length (1) /= N then
+         raise Constraint_Error with "matrices have unequal number of rows";
       end if;
 
-      --  ??? Should solve directly, is faster and more accurate
+      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;
 
-      B := Inverse (A) * X;
-   end Solve;
+         for K in MB'Range (2) loop
+            MB (MB'First (1) + J, K) := X (X'First (1) + J, K);
+         end loop;
+      end loop;
 
-   function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector is
-      B : Real_Vector (A'Range (2));
-   begin
-      Solve (A, X, B);
-      return B;
+      Forward_Eliminate (MA, MB, Det);
+      Back_Substitute (MA, MB);
+
+      return MB;
    end Solve;
 
-   function Solve (A, X : Real_Matrix) return Real_Matrix is
-      B : Real_Matrix (A'Range (2), X'Range (2));
+   ----------------------
+   -- Sort_Eigensystem --
+   ----------------------
+
+   procedure Sort_Eigensystem
+     (Values  : in out Real_Vector;
+      Vectors : in out Real_Matrix)
+   is
+
+      procedure Swap (Left, Right : Integer);
+      --  Swap Values (Left) with Values (Right), and also swap the
+      --  corresponding eigenvectors. Note that lowerbounds may differ.
+
+      function Less (Left, Right : Integer) return Boolean is
+        (Values (Left) > Values (Right));
+      --  Sort by decreasing eigenvalue, see RM G.3.1 (76).
+
+      procedure Sort is new Generic_Anonymous_Array_Sort (Integer);
+      --  Sorts eigenvalues and eigenvectors by decreasing value
+
+      procedure Swap (Left, Right : Integer) is
+      begin
+         Swap (Values (Left), Values (Right));
+         Swap_Column (Vectors, Left - Values'First + Vectors'First (2),
+                               Right - Values'First + Vectors'First (2));
+      end Swap;
+
    begin
-      Solve (A, X, B);
-      return B;
-   end Solve;
+      Sort (Values'First, Values'Last);
+   end Sort_Eigensystem;
 
    ---------------
    -- Transpose --
Index: s-gearop.adb
===================================================================
--- s-gearop.adb	(revision 178155)
+++ s-gearop.adb	(working copy)
@@ -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- --
@@ -43,6 +43,27 @@ 
       First : Integer) return Integer;
    pragma Inline_Always (Check_Unit_Last);
 
+   --------------
+   -- Diagonal --
+   --------------
+
+   function Diagonal (A : Matrix) return Vector is
+
+      N : constant Natural := Natural'Min (A'Length (1), A'Length (2));
+      R : Vector (A'First (1) .. A'First (1) + N - 1);
+
+   begin
+      for J in 0 .. N - 1 loop
+         R (R'First + J) := A (A'First (1) + J, A'First (2) + J);
+      end loop;
+
+      return R;
+   end Diagonal;
+
+   --------------------------
+   -- Square_Matrix_Length --
+   --------------------------
+
    function Square_Matrix_Length (A : Matrix) return Natural is
    begin
       if A'Length (1) /= A'Length (2) then
@@ -73,6 +94,196 @@ 
       return First + (Order - 1);
    end Check_Unit_Last;
 
+   ---------------------
+   -- Back_Substitute --
+   ---------------------
+
+   procedure Back_Substitute (M, N : in out Matrix) is
+      pragma Assert (M'First (1) = N'First (1) and then
+                     M'Last  (1) = N'Last (1));
+      Max_Col : Integer := M'Last (2);
+
+      procedure Sub_Row
+        (M      : in out Matrix;
+         Target : Integer;
+         Source : Integer;
+         Factor : Scalar);
+
+      procedure Sub_Row
+        (M      : in out Matrix;
+         Target : Integer;
+         Source : Integer;
+         Factor : Scalar) is
+      begin
+         for J in M'Range (2) loop
+            M (Target, J) := M (Target, J) - Factor * M (Source, J);
+         end loop;
+      end Sub_Row;
+
+   begin
+      for Row in reverse M'Range (1) loop
+         Find_Non_Zero : for Col in M'First (2) .. Max_Col loop
+            if Is_Non_Zero (M (Row, Col)) then
+               --  Found first non-zero element, so subtract a multiple
+               --  of this row from all higher rows, to reduce all other
+               --  elements in this column to zero.
+
+               for J in M'First (1) .. Row - 1 loop
+                  Sub_Row (N, J, Row, (M (J, Col) / M (Row, Col)));
+                  Sub_Row (M, J, Row, (M (J, Col) / M (Row, Col)));
+               end loop;
+
+               Max_Col := Col - 1;
+               exit Find_Non_Zero;
+            end if;
+         end loop Find_Non_Zero;
+      end loop;
+   end Back_Substitute;
+
+   -----------------------
+   -- Forward_Eliminate --
+   -----------------------
+
+   procedure Forward_Eliminate
+     (M   : in out Matrix;
+      N   : in out Matrix;
+      Det : out Scalar)
+   is
+      pragma Assert (M'First (1) = N'First (1) and then
+                     M'Last  (1) = N'Last (1));
+
+      function "abs" (X : Scalar) return Scalar is
+        (if X < Zero then Zero - X else X);
+
+      procedure Sub_Row
+        (M : in out Matrix;
+         Target : Integer;
+         Source : Integer;
+         Factor : Scalar);
+
+      procedure Divide_Row
+        (M, N  : in out Matrix;
+         Row   : Integer;
+         Scale : Scalar);
+
+      procedure Switch_Row
+        (M, N  : in out Matrix;
+         Row_1 : Integer;
+         Row_2 : Integer);
+
+      -------------
+      -- Sub_Row --
+      -------------
+
+      procedure Sub_Row
+        (M      : in out Matrix;
+         Target : Integer;
+         Source : Integer;
+         Factor : Scalar) is
+      begin
+         for J in M'Range (2) loop
+            M (Target, J) := M (Target, J) - Factor * M (Source, J);
+         end loop;
+      end Sub_Row;
+
+      ----------------
+      -- Divide_Row --
+      ----------------
+
+      procedure Divide_Row
+        (M, N  : in out Matrix;
+         Row   : Integer;
+         Scale : Scalar)
+      is
+      begin
+         Det := Det * Scale;
+
+         for J in M'Range (2) loop
+            M (Row, J) := M (Row, J) / Scale;
+         end loop;
+
+         for J in N'Range (2) loop
+            N (Row - M'First (1) + N'First (1), J)
+               := N (Row - M'First (1) + N'First (1), J) / Scale;
+         end loop;
+      end Divide_Row;
+
+      ----------------
+      -- Switch_Row --
+      ----------------
+
+      procedure Switch_Row
+        (M, N  : in out Matrix;
+         Row_1 : Integer;
+         Row_2 : Integer)
+      is
+         procedure Swap (X, Y : in out Scalar);
+         --  Exchange the values of X and Y
+
+         procedure Swap (X, Y : in out Scalar) is
+            T : constant Scalar := X;
+         begin
+            X := Y;
+            Y := T;
+         end Swap;
+
+      begin
+         if Row_1 /= Row_2 then
+            Det := Zero - Det;
+
+            for J in M'Range (2) loop
+               Swap (M (Row_1, J), M (Row_2, J));
+            end loop;
+
+            for J in N'Range (2) loop
+               Swap (N (Row_1 - M'First (1) + N'First (1), J),
+                     N (Row_2 - M'First (1) + N'First (1), J));
+            end loop;
+         end if;
+      end Switch_Row;
+
+      I   : Integer := M'First (1);
+
+   begin -- Forward_Eliminate
+      Det := One;
+
+      for J in M'Range (2) loop
+         declare
+            Max_I   : Integer := I;
+            Max_Abs : Scalar := Zero;
+         begin
+            --  Find best pivot in column J, starting in row I.
+            for K in I .. M'Last (1) loop
+               declare
+                  New_Abs : constant Scalar := abs M (K, J);
+               begin
+                  if Max_Abs < New_Abs then
+                     Max_Abs := New_Abs;
+                     Max_I := K;
+                  end if;
+               end;
+            end loop;
+
+            if Zero < Max_Abs then
+               Switch_Row (M, N, I, Max_I);
+               Divide_Row (M, N, I, M (I, J));
+
+               for U in I + 1 .. M'Last (1) loop
+                  Sub_Row (N, U, I, M (U, J));
+                  Sub_Row (M, U, I, M (U, J));
+               end loop;
+
+               exit when I >= M'Last (1);
+
+               I := I + 1;
+
+            else
+               Det := Zero; --  Zero, but we don't have literals
+            end if;
+         end;
+      end loop;
+   end Forward_Eliminate;
+
    -------------------
    -- Inner_Product --
    -------------------
@@ -97,6 +308,15 @@ 
       return R;
    end Inner_Product;
 
+   -------------
+   -- L2_Norm --
+   -------------
+
+   function L2_Norm (X : Vector) return Scalar is
+   begin
+      return Sqrt (Inner_Product (X, X));
+   end L2_Norm;
+
    ----------------------------------
    -- Matrix_Elementwise_Operation --
    ----------------------------------
@@ -402,6 +622,20 @@ 
       return R;
    end Outer_Product;
 
+   -----------------
+   -- Swap_Column --
+   -----------------
+
+   procedure Swap_Column (A : in out Matrix; Left, Right : Integer) is
+      Temp : Scalar;
+   begin
+      for J in A'Range (1) loop
+         Temp := A (J, Left);
+         A (J, Left) := A (J, Right);
+         A (J, Right) := Temp;
+      end loop;
+   end Swap_Column;
+
    ---------------
    -- Transpose --
    ---------------
Index: s-gearop.ads
===================================================================
--- s-gearop.ads	(revision 178155)
+++ s-gearop.ads	(working copy)
@@ -6,7 +6,7 @@ 
 --                                                                          --
 --                                 S p e c                                  --
 --                                                                          --
---          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- --
@@ -32,6 +32,50 @@ 
 package System.Generic_Array_Operations is
 pragma Pure (Generic_Array_Operations);
 
+   ---------------------
+   -- Back_Substitute --
+   ---------------------
+
+   generic
+      type Scalar is private;
+      type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+      with function "-" (Left, Right : Scalar) return Scalar is <>;
+      with function "*" (Left, Right : Scalar) return Scalar is <>;
+      with function "/" (Left, Right : Scalar) return Scalar is <>;
+      with function Is_Non_Zero (X : Scalar) return Boolean is <>;
+   procedure Back_Substitute (M, N : in out Matrix);
+
+   --------------
+   -- Diagonal --
+   --------------
+
+   generic
+      type Scalar is private;
+      type Vector is array (Integer range <>) of Scalar;
+      type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+   function Diagonal (A : Matrix) return Vector;
+
+   -----------------------
+   -- Forward_Eliminate --
+   -----------------------
+
+   --  Use elementary row operations to put square matrix M in row echolon
+   --  form. Identical row operations are performed on matrix N, must have the
+   --  same number of rows as M.
+
+   generic
+      type Scalar is private;
+      type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+      with function "-" (Left, Right : Scalar) return Scalar is <>;
+      with function "*" (Left, Right : Scalar) return Scalar is <>;
+      with function "/" (Left, Right : Scalar) return Scalar is <>;
+      with function "<" (Left, Right : Scalar) return Boolean is <>;
+      Zero, One : Scalar;
+   procedure Forward_Eliminate
+     (M   : in out Matrix;
+      N   : in out Matrix;
+      Det : out Scalar);
+
    --------------------------
    -- Square_Matrix_Length --
    --------------------------
@@ -242,6 +286,17 @@ 
      (Left  : Left_Vector;
       Right : Right_Vector) return Result_Scalar;
 
+   -------------
+   -- L2_Norm --
+   -------------
+
+   generic
+      type Scalar is private;
+      type Vector is array (Integer range <>) of Scalar;
+      with function Inner_Product (Left, Right : Vector) return Scalar is <>;
+      with function Sqrt (X : Scalar) return Scalar is <>;
+   function L2_Norm (X : Vector) return Scalar;
+
    -------------------
    -- Outer_Product --
    -------------------
@@ -332,6 +387,15 @@ 
      (Left  : Left_Matrix;
       Right : Right_Matrix) return Result_Matrix;
 
+   -----------------
+   -- Swap_Column --
+   -----------------
+
+   generic
+      type Scalar is private;
+      type Matrix is array (Integer range <>, Integer range <>) of Scalar;
+   procedure Swap_Column (A : in out Matrix; Left, Right : Integer);
+
    ---------------
    -- Transpose --
    ---------------