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 bo...@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.
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.
+
