The Computer Language
23.03 Benchmarks Game

spectral-norm Ada 2012 GNAT #3 program

source code


-- The Computer Language Benchmarks Game
-- https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
--
-- Contributed by Jim Rogers
-- Modified by Jonathan Parker (Oct 2009)

pragma Restrictions (No_Abort_Statements);
pragma Restrictions (Max_Asynchronous_Select_Nesting => 0);

with Ada.Text_Io;
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Command_Line; use Ada.Command_Line;
with Spectral_Utils;

procedure SpectralNorm is

   type Real is digits 15;

   No_of_Cores_to_Use : constant := 4;

   package Real_IO is new Ada.Text_Io.Float_Io(Real);
   package Real_Funcs is new Ada.Numerics.Generic_Elementary_Functions(Real);
   use Real_Funcs;

   N : Natural := 100;
   Vbv, Vv : Real := 0.0;
begin
   if Argument_Count = 1 then
      N := Natural'Value (Argument(1));
   end if;

   declare
      package Spectral_Utilities is new Spectral_Utils
        (Real, No_of_Tasks => No_of_Cores_to_Use, Matrix_Size => N);
      use Spectral_Utilities;

      U : Matrix := (Others => 1.0);
      V : Matrix := (Others => 0.0);
   begin
      for I in 1 .. 10 loop
         Eval_Ata_Times_U(U, V);
         Eval_Ata_Times_U(V, U);
      end loop;
      for I in V'Range loop
         Vbv := Vbv + U(I) * V(I);
         Vv  := Vv  + V(I) * V(I);
      end loop;
   end;
   Real_IO.Put(Item => Sqrt(Vbv/Vv), Fore => 1, Aft => 9, Exp => 0);
   Ada.Text_Io.New_Line;
end SpectralNorm;

generic

   type Real is digits <>;
   No_Of_Tasks : Positive;
   Matrix_Size : Positive;

package Spectral_Utils is

   type Matrix is array(Natural range 0 .. Matrix_Size-1) of Real;

   --  Evaluate matrix A at indices I, J.
 
   function Eval_A(I, J : Natural) return Real;

   --  Get   A_transpose_A_times_U = A_transpose * A * U. 
 
   procedure Eval_Ata_Times_U
     (U                     : in Matrix; 
      A_transpose_A_times_U : out Matrix);

   --  Get   AU = A * U.  Calculate only AU(Start .. Finish).
 
   procedure Eval_A_Times 
     (U      : in  Matrix; 
      Start  : in  Natural; 
      Finish : in  Natural; 
      AU     : out Matrix);

   --  Get   AU = A_transpose * U.   Calculate only AU(Start .. Finish).
 
   procedure Eval_At_Times
     (U      : in  Matrix; 
      Start  : in  Natural; 
      Finish : in  Natural; 
      AU     : out Matrix);

   pragma Inline (Eval_A_Times, Eval_At_Times);
   pragma Inline (Eval_A, Eval_Ata_Times_U);

end Spectral_Utils;

package body Spectral_Utils is

   function Eval_A (I, J : in Natural) return Real is
      Denom : constant Real := Real (((I + J) * (I + J + 1)) / 2 + I + 1);
   begin
      return 1.0 / Denom;
   end Eval_A;

   type A_Element_Pair is array (0 .. 1) of Real;

   -- Evaluate matrix A twice - at (I,J) and (I,J+1):

   function Eval_A_Twice (I, J : in Integer) return A_Element_Pair is
      Denom_0 : constant Real := Real ((I + J    )*(I + J + 1)/2 + I + 1);
      Denom_1 : constant Real := Real ((I + J + 1)*(I + J + 2)/2 + I + 1);
   begin
      return (1.0 / Denom_0, 1.0 / Denom_1);
   end Eval_A_Twice;

   -- Evaluate A_transpose (indices I and J swapped):

   function Eval_A_tr_Twice (I, J : in Integer) return A_Element_Pair is 
      Denom_0 : constant Real := Real ((I + J    )*(I + J + 1)/2 + J + 1);
      Denom_1 : constant Real := Real ((I + J + 1)*(I + J + 2)/2 + J + 2);
   begin
      return (1.0 / Denom_0, 1.0 / Denom_1);
   end Eval_A_tr_Twice;

   procedure Eval_A_Times
     (U      : in  Matrix; 
      Start  : in  Natural; 
      Finish : in  Natural; 
      Au     : out Matrix) 
   is
      Sum : Real;
      J_Index : Natural;
      A_Elements : A_Element_Pair;
   begin
      for I in Start .. Finish loop
         Sum := 0.0;
         for J in Natural range 0 .. U'Length/2 - 1 loop
            J_Index    := U'First + 2*J;
            A_Elements := Eval_A_Twice (I, J_Index);
            Sum := Sum + A_Elements(0)*U(J_Index) + A_Elements(1)*U(J_Index+1);
         end loop;
         if U'Length mod 2 = 1 then
            Sum := Sum + Eval_A(I, U'Last) * U(U'Last); -- J_Index := U'Last;
         end if;
         Au(I) := Sum;
      end loop;
   end Eval_A_Times;


   procedure Eval_At_Times
     (U      : in  Matrix; 
      Start  : in  Natural; 
      Finish : in  Natural; 
      Au     : out Matrix) 
   is
      Sum : Real;
      J_Index : Natural;
      A_Elements : A_Element_Pair;
   begin
      for I in Start .. Finish loop
         Sum := 0.0;
         for J in Natural range 0 .. U'Length/2 - 1 loop
            J_Index    := U'First + 2*J;
            A_Elements := Eval_A_tr_Twice (I, J_Index);
            Sum := Sum + A_Elements(0)*U(J_Index) + A_Elements(1)*U(J_Index+1);
         end loop;
         if U'Length mod 2 = 1 then
            Sum := Sum + Eval_A (U'Last, I) * U(U'Last); -- J_Index := U'Last;
         end if;
         Au(I) := Sum;
      end loop;
   end Eval_At_Times;


   -- Calculate A * U

   task type Matrix_A_times_U is
      pragma Storage_Size (2**20);
      entry Multiply (U : in Matrix; Start : in Natural; Finish : in Natural);
      entry Result (Start : out Natural; Finish : out Natural; R : out Matrix);
   end Matrix_A_times_U;

   task body Matrix_A_times_U is
      I1, I2 : Natural;
      AU, U_local : Matrix;
   begin
     loop
     select
        accept Multiply (U : in Matrix; Start : in Natural; Finish : in Natural) do
           I1 := Start;
           I2 := Finish;
           U_local := U;
        end Multiply;
  
        Eval_A_Times (U_local, I1, I2, AU); -- updates AU(I1..I2)
  
        accept Result (Start : out Natural; Finish : out Natural; R : out Matrix) do
           Start  := I1;
           Finish := I2;
           R(Start .. Finish) := AU(Start .. Finish);
        end Result;
     or
        terminate;
     end select;
     end loop;
   end Matrix_A_times_U;


   -- Calculate A_transpose * V

   task type Matrix_A_tr_times_V is
      pragma Storage_Size (2**20);
      entry Multiply (V : in Matrix; Start : in Natural; Finish : in Natural);
      entry Result (Start : out Natural; Finish : out Natural; R : out Matrix);
   end Matrix_A_tr_times_V;

   task body Matrix_A_tr_times_V is
      I1, I2 : Natural;
      AV, V_local : Matrix;
   begin
     loop
     select
        accept Multiply (V : in Matrix; Start : in Natural; Finish : in Natural) do
           I1 := Start;
           I2 := Finish;
           V_local := V;
        end Multiply;
  
        Eval_At_Times (V_local, I1, I2, AV);  -- AV = A_transpose * V_local
  
        accept Result (Start : out Natural; Finish : out Natural; R : out Matrix) do
           Start  := I1;
           Finish := I2;
           R(Start .. Finish) := AV(Start .. Finish);
        end Result;
     or
        terminate;
     end select;
     end loop;
   end Matrix_A_tr_times_V;


   -- Create (No_Of_Tasks-1) tasks. The final task is the environmental task, 
   -- which does its fair share of the work in procedure Eval_Ata_Times_U.

   subtype Task_Range is Positive range 1 .. No_Of_Tasks-1;

   Partial_Matrix_A_times_U    : array (Task_Range) of Matrix_A_times_U; 
   Partial_Matrix_A_tr_times_V : array (Task_Range) of Matrix_A_tr_times_V; 


   procedure Eval_Ata_Times_U 
     (U                     : in  Matrix; 
      A_transpose_A_times_U : out Matrix) 
   is
      V, Partial_Product : Matrix;
   
      Segment_Length : constant Integer := U'Length / No_Of_Tasks + 1;
      -- Gives the 1st few tasks a slightly greater share of the work.
 
      I1, I2, J1, J2 : Natural;
   begin
      I1 := V'First;
      I2 := V'First + Segment_Length - 1;
      I2 := Integer'Min (I2, V'Last); 

      -- Start running the tasks in Task_Range:

      for k in Task_Range loop
         Partial_Matrix_A_times_U(k).Multiply (U, I1, I2);
         I1 := I2 + 1;
         I2 := I2 + Segment_Length;
         I2 := Integer'Min (I2, V'Last);
      end loop;

      Eval_A_Times (U, I1, V'Last, V); -- Env task updates V(I1 .. V'Last).

      -- Rendezvous with tasks to get partial results. Write results to V:

      for k in Task_Range loop
         Partial_Matrix_A_times_U(k).Result (J1, J2, Partial_Product);
         V(J1 .. J2) := Partial_Product(J1 .. J2);
      end loop;

      -- The result, stored in V, is A*U. Next get A_transpose * (A*U).

      I1 := V'First;
      I2 := V'First + Segment_Length - 1;
      I2 := Integer'Min (I2, V'Last); 

      for k in Task_Range loop
         Partial_Matrix_A_tr_times_V(k).Multiply (V, I1, I2);
         I1 := I2 + 1;
         I2 := I2 + Segment_Length;
         I2 := Integer'Min (I2, V'Last);
      end loop;

      Eval_At_Times (V, I1, V'Last, A_transpose_A_times_U); 
      -- Env. task updates A_transpose_A_times_U (I1 .. V'Last).

      for k in Task_Range loop
         Partial_Matrix_A_tr_times_V(k).Result (J1, J2, Partial_Product);
         A_transpose_A_times_U(J1 .. J2) := Partial_Product(J1 .. J2);
      end loop;

   end Eval_Ata_Times_U;

end Spectral_Utils;


    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
GNATMAKE 11.3.0


Wed, 25 Jan 2023 03:32:34 GMT

MAKE:
gnatchop -r -w spectralnorm.gnat-3.gnat
splitting spectralnorm.gnat-3.gnat into:
   spectralnorm.adb
   spectral_utils.ads
   spectral_utils.adb
gnatmake -O3 -fomit-frame-pointer -march=ivybridge -gnatNp -f spectralnorm.adb -o spectralnorm.gnat-3.gnat_run 
x86_64-linux-gnu-gcc-11 -c -O3 -fomit-frame-pointer -march=ivybridge -gnatNp spectralnorm.adb
x86_64-linux-gnu-gcc-11 -c -O3 -fomit-frame-pointer -march=ivybridge -gnatNp spectral_utils.adb
x86_64-linux-gnu-gnatbind-11 -x spectralnorm.ali
x86_64-linux-gnu-gnatlink-11 spectralnorm.ali -O3 -fomit-frame-pointer -march=ivybridge -o spectralnorm.gnat-3.gnat_run

3.23s to complete and log all make actions

COMMAND LINE:
./spectralnorm.gnat-3.gnat_run 5500

PROGRAM OUTPUT:
1.274224153