Re: Interfacing Ada multidimensional arrays with Fortran.

[Simon Wright <simon@pushface.org>]

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).