A thicker binding for Lapack

A thicker binding for Lapack

[jpwoodruff]

There is a complete thin binding for lapack available which was
evidently started by Wasu Chaopanon in about 1996 and was recently
finished by Nasser Abbasi.  There about 1000 subroutines; their
parameter signatures exactly translate the Fortran, and each procedure
is Pragma Imported.

I've been considering adding value to that binding to address two
perceived problems:

The specification involving the integer index used for matrix types is
inconsistent between the binding and Ada.Numerics.Generic_Real_Arrays.
The binding declares that type because the user's code must comply
with Fortran's integer.  The Ada compiler is obliged to use that same
type.  Unfortunately Generic_Real_Arrays uses the base type integer
for index, so there is no expectation that the representation would be
the same as Fortran-compiled library.

Another issue to resolve is that lapack, being Fortran, defines
column-major storage while Ada is row-major.  Parameters passing
between these two conventions evidently must be transposed.

I think I have a technique for interoperating with the Ada.Numerics
types.  My plan involves "wrapping" the lapack routines in a layer
that converts types and transposes matrices.

Here is a typical procedure specification as translated directly from
the fortran, with my transformation to Ada standard types.

   procedure SGESV
     (N    : Natural;
      NRHS : Natural;
      A    : in out Real_Arrays.Real_Matrix;   -- instance of Generic_
      LDA  : Natural ;
      IPIV : out Integer_Vector;
      B    : in out Real_Arrays.Real_Matrix;
      LDB  : Natural ;
      INFO : out Integer) ;

The math requirements for this procedure are very much like the
"Solve" that was discussed in a recent thread.

   -- function Solve (A : Real_Matrix; X : Real_Vector) return
   -- Real_Vector;

The B parameter in sgesv resembles X in the "Why X as argument name?"
discussion.  I am trying to make the lapack item as serviceable as the
gnat standard item.

My trials at maintaining parameter definitions that echo the Fortran
binding have led only to grief.  The integer inputs N, NRHS, LDA, LDB
make precious little sense, considering the transposition of the
matrices that takes place inside the body of procedure sgesv.

After several drafts I'm thinking that a thick binding might look
like this:

   procedure SGESV
      (A   : in out Real_Arrays.Real_Matrix;
      IPIV : out Integer_Vector;
      B    : in out Real_Arrays.Real_Matrix;
      INFO : out Integer) ;

Programming with Ada attributes seems to clarify the relations between
the several parameters.  Now I think that the integer parameters will
be computed inside the body of the overlay procedure.  Still working
on the details.

Guided by the spiral model I'll try a couple examples.  Later if there
is a benefit to Ada programmers interested in numeric analysis I'll
address the scale of the massive lapack library.

Suggestions are welcome!

John

Re: A thicker binding for Lapack

[Shark8]

On Thursday, December 27, 2012 12:48:14 PM UTC-6, jpwoo...@gmail.com wrote:
>
> Another issue to resolve is that lapack, being Fortran, defines
> column-major storage while Ada is row-major.  Parameters passing
> between these two conventions evidently must be transposed.

This is a non-issue; Ada allows for the specification of row- or column-major ordering via the Convention pragma (or aspect).

Example:
   B : Array( Integer Range <>, Integer Range <> ) of Integer:=
     (1..5 => (1..5 => 0) );
   Pragma Convention( Fortran, B ); -- Forces B to use col-major order.

-- NOTE: this pragma (or aspect) can be attached to the TYPE definition.
   Type Test_Array is Array( Integer Range <>, Integer Range <> ) of Integer
        with	Convention => Fortran;
   
   C : Test_Array:= (1..5 => (1..5 => 0) );

Re: A thicker binding for Lapack

[Simon Wright]

jpwoodruff@gmail.com writes:

> My trials at maintaining parameter definitions that echo the Fortran
> binding have led only to grief.  The integer inputs N, NRHS, LDA, LDB
> make precious little sense, considering the transposition of the
> matrices that takes place inside the body of procedure sgesv.

Completely agree with this.

> After several drafts I'm thinking that a thick binding might look
> like this:
>
>    procedure SGESV
>       (A   : in out Real_Arrays.Real_Matrix;
>       IPIV : out Integer_Vector;
>       B    : in out Real_Arrays.Real_Matrix;
>       INFO : out Integer) ;

[see extract from SGESV man page below]

There may be some deep reason for making A in-out, but if you don't see
a need to preserve the internal LU decomposition I'd make a copy
internally, and also keep IPIV internal.

Also, I'd make B 'in' and use a third 'out' parameter Result.

And, perhaps INFO should be replaced by Constraint_Error?

[from SGESV man page]

PURPOSE
       SGESV computes the solution to a real system of linear equations
          A * X = B, where A is an N-by-N matrix and X  and  B  are  N-by-NRHS
       matrices.   The  LU  decomposition with partial pivoting and row inter-
       changes is used to factor A as
          A = P * L * U,
       where P is a permutation matrix, L is unit lower triangular, and  U  is
       upper  triangular.   The  factored  form of A is then used to solve the
       system of equations A * X = B.

ARGUMENTS
       [...]

       A       (input/output) REAL array, dimension (LDA,N)
               On entry, the N-by-N coefficient matrix A.  On exit,  the  fac-
               tors  L and U from the factorization A = P*L*U; the unit diago-
               nal elements of L are not stored.

       [...]

       IPIV    (output) INTEGER array, dimension (N)
               The pivot indices that define the permutation matrix P;  row  i
               of the matrix was interchanged with row IPIV(i).

       B       (input/output) REAL array, dimension (LDB,NRHS)
               On entry, the N-by-NRHS matrix of right hand side matrix B.  On
               exit, if INFO = 0, the N-by-NRHS solution matrix X.

       [...]

       INFO    (output) INTEGER
               = 0:  successful exit
               < 0:  if INFO = -i, the i-th argument had an illegal value
               > 0:  if INFO = i, U(i,i) is exactly zero.   The  factorization
               has  been  completed,  but the factor U is exactly singular, so
               the solution could not be computed.

    

Re: A thicker binding for Lapack

[jpwoodruff]

On Thursday, December 27, 2012 12:18:09 PM UTC-7, Shark8 wrote:
> On Thursday, December 27, 2012 12:48:14 PM UTC-6, jpwoo...@gmail.com wrote:
> 
> >
> 
> > Another issue to resolve is that lapack, being Fortran, defines
> 
> > column-major storage while Ada is row-major.  Parameters passing
> 
> > between these two conventions evidently must be transposed.
> 
> 
> 
> This is a non-issue; Ada allows for the specification of row- or column-major ordering via the Convention pragma (or aspect).
> 
The use case I'm working on is to allow a mathematical program that
already uses matrices to be enhanced with lapack operators on the same
type of operand.  

If I'm not mistaken, pragma convention will affect matrices that are
used on the lapack interface, but other matrix types in the program
remain row-wise.  

This is the inconsistency I'm working to avoid unless every matrix
type is conventioned, and it's not clear how to catch them all.

The alignment is visible to the user through (at least) the program's
I/O.

Conclusion:   transpose appears necessary.

John

Re: A thicker binding for Lapack

[Simon Wright]

jpwoodruff@gmail.com writes:

> Conclusion:   transpose appears necessary.

For info, the compiler (well, GNAT) will automatically transpose when
doing assignment between arrays with different conventions. And it's
quicker (not by very much at -O2, though).

with Ada.Calendar; use Ada.Calendar;
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Float_Random;
with Interfaces.Fortran;
with System.Generic_Array_Operations;
procedure Sauvage_Timing is

   package BLAS is

      --  Copied from old GNAT Interfaces.Fortran.BLAS.

      --  Vector types

      type Real_Vector is array (Integer range <>)
        of Interfaces.Fortran.Real;

      type Complex_Vector is array (Integer range <>)
        of Interfaces.Fortran.Complex;

      type Double_Precision_Vector is array (Integer range <>)
        of Interfaces.Fortran.Double_Precision;

      type Double_Complex_Vector is array (Integer range <>)
        of Interfaces.Fortran.Double_Complex;

      --  Matrix types

      type Real_Matrix is array (Integer range <>, Integer range <>)
        of Interfaces.Fortran.Real;

      type Double_Precision_Matrix
         is array (Integer range <>, Integer range <>)
        of Interfaces.Fortran.Double_Precision;

      type Complex_Matrix is array (Integer range <>, Integer range <>)
        of Interfaces.Fortran.Complex;

      type Double_Complex_Matrix is array (Integer range <>, Integer range <>)
        of Interfaces.Fortran.Double_Complex;

   end BLAS;

   procedure Transpose
   is new System.Generic_Array_Operations.Transpose
     (Scalar => Interfaces.Fortran.Double_Precision'Base,
      Matrix => BLAS.Double_Precision_Matrix);

   type Double_Precision_Matrix is array (Integer range <>, Integer range <>)
     of Interfaces.Fortran.Double_Precision;
   pragma Convention (Fortran, Double_Precision_Matrix);

   A, B : BLAS.Double_Precision_Matrix (1 .. 100, 1 .. 100);
   C : Double_Precision_Matrix (1 .. 100, 1 .. 100);
   pragma Volatile (B);
   pragma Volatile (C);

   Gen : Ada.Numerics.Float_Random.Generator;

   Start, Finish : Time;

   use type Interfaces.Fortran.Double_Precision;

begin

   Ada.Numerics.Float_Random.Reset (Gen);

   for J in A'Range (1) loop
      for K in A'Range (2) loop
         A (J, K) :=
           Interfaces.Fortran.Double_Precision (J * K)
           * Interfaces.Fortran.Double_Precision
           (Ada.Numerics.Float_Random.Random (Gen));
      end loop;
   end loop;

   Start := Clock;
   for J in 1 .. 100 loop
      Transpose (A, B);
   end loop;
   Finish := Clock;
   Put_Line ("Transpose took " & Duration'Image (Finish - Start));

   Start := Clock;
   for J in 1 .. 100 loop
      C := Double_Precision_Matrix (A);
   end loop;
   Finish := Clock;
   Put_Line ("Assignment took" & Duration'Image (Finish - Start));

   declare
      Same : Boolean := True;
   begin
      for J in 1 .. 100 loop
         for K in 1 .. 100 loop
            if B (J, K) /= C (K, J) then
               Same := False;
            end if;
         end loop;
      end loop;
      Put_Line ("B = ~C: " & Same'Img);
   end;

end Sauvage_Timing;

Re: A thicker binding for Lapack

[jpwoodruff]

On Saturday, December 29, 2012 12:59:46 PM UTC-7, Simon Wright wrote:
> jpwoodruff@gmail.com writes:
> 
> 
> 
> > Conclusion:   transpose appears necessary.
> 
> 
> 
> For info, the compiler (well, GNAT) will automatically transpose when
> 
> doing assignment between arrays with different conventions. And it's
> 
> quicker (not by very much at -O2, though).
> 

I'll be darned.  You are both right about pragma Convention.  I am
going to unwrite some code.  Even my old gnat gpl 2010 got it right. 
I need to (re)learn to trust Lady Ada.


From Simon's context, I found System.Generic_Real_LAPACK. There is
quite a lot there.  I want to spend some time reading that, then
determine if (as I suspect) the problem I stated has already been
solved by the structures there.

John

Re: A thicker binding for Lapack

[Simon Wright]

jpwoodruff@gmail.com writes:

> On Saturday, December 29, 2012 12:59:46 PM UTC-7, Simon Wright wrote:
>> jpwoodruff@gmail.com writes:
>> 
>> > Conclusion:   transpose appears necessary.
>> 
>> For info, the compiler (well, GNAT) will automatically transpose when
>> doing assignment between arrays with different conventions. And it's
>> quicker (not by very much at -O2, though).
>> 
>
> I'll be darned.  You are both right about pragma Convention.  I am
> going to unwrite some code.  Even my old gnat gpl 2010 got it right. 
> I need to (re)learn to trust Lady Ada.
>
>
> From Simon's context, I found System.Generic_Real_LAPACK. There is
> quite a lot there.  I want to spend some time reading that, then
> determine if (as I suspect) the problem I stated has already been
> solved by the structures there.

System.Generic_Real_LAPACK looks a bit off to me - lots of unchecked
conversions of address-of-array into access-to-another-array-type. But
it may do what you need.

NB, GNATs later than GCC 4.6 or GNAT GPL 2010 don't use LAPACK or
BLAS. I wrote a bit about this on comp.lang.ada at [1].

[1] http://coding.derkeiler.com/Archive/Ada/comp.lang.ada/2012-07/msg00138.html

Re: A thicker binding for Lapack

[jpwoodruff]

On Sunday, December 30, 2012 12:41:54 PM UTC-7, Simon Wright wrote:
> jpwoodruff@gmail.com writes:
> 
> 
> 
> > On Saturday, December 29, 2012 12:59:46 PM UTC-7, Simon Wright wrote:
> 
> >> jpwoodruff@gmail.com writes:
> 
> >> 
> 
> >> > Conclusion:   transpose appears necessary.
> 
> >> 
> 
> >> For info, the compiler (well, GNAT) will automatically transpose when
> 
> >> doing assignment between arrays with different conventions. And it's
> 
> >> quicker (not by very much at -O2, though).
> 
> >> 
> 
> >
> 
> > I'll be darned.  You are both right about pragma Convention.  I am
> 
> > going to unwrite some code.  Even my old gnat gpl 2010 got it right. 
> 
> > I need to (re)learn to trust Lady Ada.
> 
> >
> 
> >
> 
> > From Simon's context, I found System.Generic_Real_LAPACK. There is
> 
> > quite a lot there.  I want to spend some time reading that, then
> 
> > determine if (as I suspect) the problem I stated has already been
> 
> > solved by the structures there.
> 
> 
> 
> System.Generic_Real_LAPACK looks a bit off to me - lots of unchecked
> 
> conversions of address-of-array into access-to-another-array-type. But
> 
> it may do what you need.
> 
> 
> 
> NB, GNATs later than GCC 4.6 or GNAT GPL 2010 don't use LAPACK or
> 
> BLAS. I wrote a bit about this on comp.lang.ada at [1].
> 
> 
> 
> [1] http://coding.derkeiler.com/Archive/Ada/comp.lang.ada/2012-07/msg00138.html

I'm uncertain about what Ada lapack materials are in play.  Please
help me get the story of the several variations straight.

about System.Generic_Real_LAPACK -- 
Is it true that this structure that I see in my gpl 2010 is obsolete?
Has it been removed in later gnats?  Has interfaces.fortran.lapack met
the same fate?

About lapack.ads on sourceforge, attributed to NMA and SJW: --
Is this in fact the only low-level binding to the Fortran linear
algebra subroutines?

Properties of lapack.ads:
 - largely auto-generated.
 - covers the entire lapack library: single, double precision, real
and complex.  
 - defines array types with components and indices specified from
types defined by interfaces.fortran.
 - uses pragma fortran for vector and matrix array types.

Evidently conversion between columnwise and row-wise storage is
handled transparently by Ada language specification of pragma.
Transparent transposition happens on assignment (but not on parameter
elaboration?).  Some testing is in order here.

I considered enumerating the properties of System.Generic_Real_LAPACK
but if it's leaving, why bother?

John