Interfacing Ada multidimensional arrays with Fortran.

Interfacing Ada multidimensional arrays with Fortran.

[David Sauvage]

Hi,

Concerning multidimensional arrays, Ada use row-major order [0] while
Fortran use column-major order [1].

It seems to exist two possibilities to handle this issue when calling
Fortran procedures from Ada, and three possibilities when calling BLAS
[2] or LAPACK [3] Fortran specific routines.

1 - Use the pragma Convention (Fortran, T_Real_Matrix)
on the array type, so that the array will be stored in column-major
order, as needed by Fortran.
Any feedbacks concerning the successful use of this pragma Convention
(Fortran, type) & GNATGPL 2010 or FSF GCC would be appreciate.

2 - Transpose the array in the Ada side when calling the Fortran code,
so that the Fortran code will load the array properly.

3 - Concerning BLAS & LAPACK routines, there are parameters to
indicates the form of the array (transposed or not) to the Fortran
code, so that the Fortran code load the array properly.

In Interfaces.Fortran.BLAS, there are Multidimensional arrays types
that are already defined ;
(Real_Matrix, Double_Precision_Matrix, Complex_Matrix,
Double_Complex_Matrix)
But the corresponding pragma Convention (Fortran, type) are not
defined for them. So any user that would re-use those types can not
use pragma Convention, and use possibility 2 or 3 above,

It would be interesting to know the story behind the scene of why
array types declared in  Interfaces.Fortran.BLAS do not use pragma
Convention (Fortran, *) ?

What are the performance issues by using the three possibilities
above ? (may be the community already had some interesting information
about this)

Some would think the pragma Convention (Fortran, type) should be more
efficient, but may be it is not that simple, as apparently the choice
has been made to avoid [4] pragma Convention (Fortran, type).


Cheers


[0] http://en.wikipedia.org/wiki/Row-major_order
[1] http://www.adaic.org/resources/add_content/standards/05rm/html/RM-B-5.html
[2] http://www.netlib.org/blas/
[3] http://www.netlib.org/lapack/
[4] i-forbla.ads

Re: Interfacing Ada multidimensional arrays with Fortran.

[Simon Wright]

David Sauvage <sauvage.david@gmail.com> writes:

> Concerning multidimensional arrays, Ada use row-major order [0] while
> Fortran use column-major order [1].

> 3 - Concerning BLAS & LAPACK routines, there are parameters to
> indicates the form of the array (transposed or not) to the Fortran
> code, so that the Fortran code load the array properly.

This is true of BLAS, but I don't believe it's true of LAPACK.

> In Interfaces.Fortran.BLAS, there are Multidimensional arrays types
> that are already defined ;
> (Real_Matrix, Double_Precision_Matrix, Complex_Matrix,
> Double_Complex_Matrix)
> But the corresponding pragma Convention (Fortran, type) are not
> defined for them. So any user that would re-use those types can not
> use pragma Convention, and use possibility 2 or 3 above,

Interfaces.Fortran.BLAS is an internal GNAT unit, ie part of GNAT's
implementation of the standard Ada.Numerics.Generic_*_Arrays, not part
of the standard itself.

> It would be interesting to know the story behind the scene of why
> array types declared in  Interfaces.Fortran.BLAS do not use pragma
> Convention (Fortran, *) ?

There I can't help you.

The implementor of Ada.Numerics.Generic_Real_Arrays has decided to
declare Interfaces.Fortran.BLAS using Ada order, so that to convert from
the Ada order required in the standard for these units he uses the
Transpose operation (note, a copy is required anyway because LAPACK
doesn't preserve the input matrix).

Functionally, he could equally well have declared in Fortran order, as
you suggest:

   with Ada.Text_IO; use Ada.Text_IO;
   with Interfaces.Fortran.BLAS;
   procedure Sauvage is
      type Real_Matrix is array (Integer range <>, Integer range <>)
        of Interfaces.Fortran.Real;
      pragma Convention (Fortran, Real_Matrix);
      Theirs : constant Interfaces.Fortran.BLAS.Real_Matrix :=
        (1 => (1 => 1.0, 2 => 2.0),
         2 => (1 => 3.0, 2 => 4.0));
      Mine : Real_Matrix (1 .. 2, 1 .. 2);
   begin
      Mine := Real_Matrix (Theirs);      -- transposition occurs here
      for J in Mine'Range (1) loop
         for K in Mine'Range (2) loop
            Put_Line (Mine (J, K)'Img);
         end loop;
      end loop;
   end Sauvage;

and it might have been quicker.

Quite how the implementor of Ada.Numerics.Generic_Complex_Arrays managed
the transposition is less than clear to me!

> What are the performance issues by using the three possibilities
> above ? (may be the community already had some interesting information
> about this)

I haven't measured this. Anyone else?

> Some would think the pragma Convention (Fortran, type) should be more
> efficient, but may be it is not that simple, as apparently the choice
> has been made to avoid [4] pragma Convention (Fortran, type).

Re: Interfacing Ada multidimensional arrays with Fortran.

[Simon Wright]

Simon Wright <simon@pushface.org> writes:

>> What are the performance issues by using the three possibilities
>> above ? (may be the community already had some interesting information
>> about this)
>
> I haven't measured this. Anyone else?

On Mac OS X, using GCC 4.6.0, with -O1, applying
System.Generic_Array_Operations.Transpose 10 times to a 100x100 matrix
took 900 us, whereas using assignment with convention Fortran took 90
us.

With -O2, the assignment was hoisted out of the loop.

I then made the target matrices volatile; the results are now

-O1:
Transposition: 1200 us
Assignment:     300 us

-O2:
Transposition:  500 us
Assignment:     300 us

Re: Interfacing Ada multidimensional arrays with Fortran.

[David Sauvage]

On May 28, 8:45 pm, Simon Wright <si...@pushface.org> wrote:
...
> -O1:
> Transposition: 1200 us
> Assignment:     300 us
>
> -O2:
> Transposition:  500 us
> Assignment:     300 us

Using this testcase [1], here are my results using Intel Atom CPU N270
@ 1.60GHz / Linux (launched with root privilege) :

 -O1:
 Transposition: 2202 us
 Assignment:    1014 us

 -O2:
 Transposition: 2556 us
 Assignment:     885 us

Transposition (using assignment via pragma Convention) seems to be
quicker than Transposition (without pragma Convention).
The reason why pragma Convention (Fortran, Type) is not used in
Interfaces.Fortran... is unknown and seems to give slower compute
time.

[1]
-- gnatmake -f compare.adb -cargs -gnat05 -O2
-- gnatmake -f compare.adb -cargs -gnat05 -O1
with Interfaces.Fortran.BLAS;

with Ada.Text_IO,
     Ada.Calendar,
     Ada.Numerics.Generic_Real_Arrays;

procedure Compare is
   Start, Stop : Ada.Calendar.Time;
   use type Ada.Calendar.Time;

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

   package GRA is new Ada.Numerics.Generic_Real_Arrays (
      Interfaces.Fortran.Real);

   Row, Column : constant Positive := 100;
   Iteration   : constant Positive := 10;

   M : Real_Matrix (1 .. Row, 1 .. Column)      := (others => (others
=> 2.0));
   pragma Volatile (M);

   MFA, MFB : GRA.Real_Matrix (1 .. Row, 1 .. Column) := (others =>
(others => 2.0));
   pragma Volatile (MFA);
   pragma Volatile (MFB);

   use type Interfaces.Fortran.Real;
begin

   Start := Ada.Calendar.Clock;
   for I in 1 .. Iteration loop
      M := Real_Matrix (MFB);
   end loop;
   Stop := Ada.Calendar.Clock;
   Ada.Text_IO.Put_Line
     ("Assignation (Transposition via pragma Convention)" &
      Duration'Image (Stop - Start));

   Start := Ada.Calendar.Clock;
   for I in 1 .. Iteration loop
      MFA := GRA.Transpose (MFB);
   end loop;
   Stop := Ada.Calendar.Clock;
   Ada.Text_IO.Put_Line
     ("Transposition" & Duration'Image (Stop - Start));

end Compare;