Tasking for Mandelbrot program

[Georg Bauhaus <rm.tsoh.plus-bug.bauhaus@maps.futureapps.de>]

A Mandelbrot program has been submitted to the language
Shootout by Jim Rogers, Pascal Obry, and
Gautier de Montmollin.  Given the no tricks approach
that it shares with other entries, it is performing
well, I'd think.  But, it is sequential.

The patch below adds tasking.

A few notes on patch:
The computation is split into parts, one per task;
many entries seem to use a similar approach.
The number of tasks is set high, task switching does
not seem to matter much in this program.
The alternative of making the environment task perform
the same subprogram as the tasks concurrently
and matching the number of tasks to the number of cores
was not significantly faster.  The alternative
spoils the simplicity of the main block, so I set
a higher number of tasks and got basically the same
performance.

Maybe these are possible improvements:

- Ada.Streams.Stream_IO? (In total absence of Text_IO,
  GNAT.IO performs well.)

- adjust compilation options to match the FPT hardware
  (type Real is now digits 16 because of this, cf.
  the spectral-norm entry at the Shootout)

Comments?  Can I ask what the authors think?


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

with Ada.Command_Line; use Ada.Command_Line;
with Ada.Numerics.Generic_Complex_Types;

with Interfaces;       use Interfaces;
with GNAT.IO;          use GNAT.IO;


procedure Mandelbrot is
   type Real is digits 16;
   package R_Complex is new Ada.Numerics.Generic_Complex_Types (Real);
   use R_Complex;
   Iter           : constant := 50;
   Limit          : constant := 4.0;
   Size           : Positive;
   Zero           : constant Complex := (0.0, 0.0);
   Two_on_Size    : Real;

   subtype Output_Queue is String;
   type Output is access Output_Queue;

   task type X_Step is
      entry Compute_Z (Y1, Y2 : Natural);
      entry Get_Output (Result : out Output; Last : out Natural);
   end X_Step;

   procedure Allocate_Output_Queue (Y1, Y2 : Natural; Result : out Output);

   procedure Compute
     (Y1, Y2 : Natural; Result : Output; Last : out Natural)
   is
      Bit_Num     : Natural    := 0;
      Byte_Acc    : Unsigned_8 := 0;
      Z, C        : Complex;
   begin
      Last := 0;
      for Y in Y1 .. Y2 - 1 loop
         for X in 0 .. Size - 1 loop
            Z := Zero;
            C := Two_on_Size * (Real (X), Real (Y)) - (1.5, 1.0);

            declare
               Z2 : Complex;
            begin
               for I in 1 .. Iter + 1 loop
                  Z2 := (Z.re ** 2, Z.im ** 2);
                  Z  := (Z2.re - Z2.im, 2.0 * Z.re * Z.im) + C;
                  exit when Z2.re + Z2.im > Limit;
               end loop;

               if Z2.re + Z2.im > Limit then
                  Byte_Acc := Shift_Left (Byte_Acc, 1) or 16#00#;
               else
                  Byte_Acc := Shift_Left (Byte_Acc, 1) or 16#01#;
               end if;
            end;

            Bit_Num := Bit_Num + 1;

            if Bit_Num = 8 then
               Last := Last + 1;
               Result (Last) := Character'Val (Byte_Acc);
               Byte_Acc := 0;
               Bit_Num  := 0;
            elsif X = Size - 1 then
               Byte_Acc := Shift_Left (Byte_Acc, 8 - (Size mod 8));
               Last := Last + 1;
               Result (Last) := Character'Val (Byte_Acc);
               Byte_Acc := 0;
               Bit_Num  := 0;
            end if;
         end loop;
      end loop;
   end Compute;

   task body X_Step is
      Data        : Output;
      Pos         : Natural;
      Y1, Y2      : Natural;
   begin
      accept Compute_Z (Y1, Y2 : Natural) do
         X_Step.Y1 := Y1;
         X_Step.Y2 := Y2;
      end Compute_Z;

      Allocate_Output_Queue (Y1, Y2, Data);
      Compute (Y1, Y2, Data, Pos);

      accept Get_Output (Result : out Output; Last : out Natural) do
         Result := Data;
         Last := Pos;
      end Get_Output;
   end X_Step;

   procedure Allocate_Output_Queue (Y1, Y2 : Natural; Result : out
Output) is
   begin
      Result := new Output_Queue (1 .. (Y2 - Y1 + 8) * Size / 8);
   end Allocate_Output_Queue;


begin
   Size := Positive'Value (Argument (1));
   Two_on_Size := 2.0 / Real (Size);

   Put_Line ("P4");
   Put_Line (Argument (1) & " " & Argument (1));

   declare
      No_Of_Workers : constant := 16;
      Chunk_Size    : constant Positive :=
        (Size + No_Of_Workers) / No_Of_Workers;
      Pool          : array (0 .. No_Of_Workers) of X_Step;
      pragma          Assert (Pool'Length * Chunk_Size >= Size);
      Buffer        : Output;
      Last          : Natural;
   begin
      pragma Assert (Pool'First = 0);

      for P in Pool'Range loop
         Pool (P).Compute_Z
           (Y1 => P * Chunk_Size,
            Y2 => Positive'Min ((P + 1) * Chunk_Size, Size));
      end loop;

      for P in Pool'Range loop
         Pool (P).Get_Output (Buffer, Last);
         Put (Buffer (Buffer'First .. Last));
      end loop;
   end;

end Mandelbrot;