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

[Simon Wright]

Post by David Sauvage
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.

Post by David Sauvage
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.

Post by David Sauvage
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!

Post by David Sauvage
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?

Post by David Sauvage
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).

[Simon Wright]

Post by Simon Wright

Post by David Sauvage
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

[David Sauvage]

On May 28, 8:45 pm, Simon Wright <***@pushface.org> wrote:
...

Post by Simon Wright
Transposition: 1200 us
Assignment:     300 us
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;