Re: Interfacing Ada to C for linear algebra

[Simon Wright <simon@pushface.org>]

"vincent.diemunsch@gmail.com" <vincent.diemunsch@gmail.com> writes:

> Hello everybody,
>
> I would like to use in Ada some linear algebra libraries written in C,
> especially the SuiteSparse :
>  http://www.cise.ufl.edu/research/sparse/SuiteSparse/
> These libraries handle Sparse Matrices that are in the format used
> internally by MatLab : a column-compressed version of the standard
> matrices.
>
> Ada defines Generic_Real_Arrays to handle vectors and matrices. So I
> would like to simply create a Generic_Real_Sparse_Arrays package to
> define MatLab compatible matrices. Since the internal format of MatLab
> matrices is in C, in need to do a conversion between my Sparse_Matrix
> and the C arrays that represent the Sparse Matrix in MatLab.
>
> Looking at the implementation of Generic_Real_Arrays in GNAT, I see
> that they internally call the BLAS and LAPACK packages, written in
> Fortran, through the use of two packages : System.Generic_Real_BLAS
> and System.Generic_Real_LAPACK. These package do exactly the same as
> what I want to do : they convert Ada arrays of Ada floating point type
> in Fortran arrays of Fortran Float type and back, without using
> copying, so that the compiler produces no code for this conversion
> (using pragma Inline_Always and so on). This is important because if
> we use optimized versions of linear algebra packages it is because we
> manipulate huge arrays and we don't want to loose time in copying
> values from on representation to another !

The LAPACK interface used in Generic_Real_Arrays deals with matrices in
Ada (row-major) format and - usually, anyway - transposes in order to
call the Fortran implementations. Unlike BLAS, LAPACK itself doesn't
have options to deal with row-major matrices.

We did some timings and found that it's quite a bit quicker to declare
the BLAS/LAPACK interfaces in terms of Fortran (column-major) matrices
(declared using pragma Convention (Fortran)) and have the compiler
implement the necessary moves on assignment, rather than using a library
Transpose subprogram. But make no mistake, the Generic_Real_Arrays code
definitely copies! (look at the body for confirmation).

> I suppose that since C is very common, there must already exist a very
> good implementation of conversion of Ada arrays of floating point
> types in C arrays of C.Float, idealy using the GNAT
> implementation. This could save me a lot of time doing interfaces...

GNAT and GCC share the same compiler backend, so the underlying numeric
types are the same. Just tell the compiler what to import:

   procedure dgeev
     (Jobv_L :        Character;
      Jobv_R :        Character;
      N      :        Positive;
      A      : in out Real_Arrays.Real_Matrix;
      Ld_A   :        Positive;
      W_R    :    out Real_Arrays.Real_Vector;
      W_I    :    out Real_Arrays.Real_Vector;
      V_L    :    out Real_Arrays.Real_Matrix;
      Ld_V_L :        Integer;
      V_R    :    out Real_Arrays.Real_Matrix;
      Ld_V_R :        Integer;
      Work   :    out Real_Arrays.Real_Vector;
      L_Work :        Integer;
      Info   :    out Integer);
   pragma Import (Fortran, dgeev, "dgeev_");

(Real_Arrays.Real_Matrix etc from Generic_Real_Arrays). The Fortran
convention means that the Character and Integer arguments are passed by
reference without your needing to fake up accesses.

More of this at https://sourceforge.net/projects/gnat-math-extn/ .