Re: A thicker binding for Lapack

[Simon Wright <simon@pushface.org>]

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;