diff --git a/.gitignore b/.gitignore index 4dfbbc405..e08b22432 100644 --- a/.gitignore +++ b/.gitignore @@ -81,10 +81,9 @@ func.c *.log3 ADA/*.ali -ADA/args -ADA/hello ADA/nstream -ADA/nstream_array +ADA/p2p +ADA/stencil ADA/transpose C1z/nstream diff --git a/ADA/CMakeLists.txt b/ADA/CMakeLists.txt new file mode 100644 index 000000000..def3f7420 --- /dev/null +++ b/ADA/CMakeLists.txt @@ -0,0 +1,67 @@ +# +# CMake port of ADA/Makefile. CMake has no native Ada language support +# (no LANGUAGES Ada for project()/enable_language()), so this drives +# gnatmake directly via add_custom_command(), the same approach used +# elsewhere in this repo for build steps CMake can't natively express +# (e.g. common/Stencil/loop_gen code generation in the C-family ports). +# +# Standalone project: `cd ADA && cmake -S . -B build` works on its own. +# Also add_subdirectory()'d from the repo root. +# + +cmake_minimum_required(VERSION 3.18) +project(PRK_ADA NONE) +enable_testing() + +find_program(GNATMAKE_EXECUTABLE gnatmake) + +if(NOT GNATMAKE_EXECUTABLE) + message(STATUS "ADA: skipped, no gnatmake (GNAT Ada compiler) found") + return() +endif() + +message(STATUS "ADA: gnatmake found at ${GNATMAKE_EXECUTABLE}") + +# prk_ada_executable() +# Compiles .adb with gnatmake into ${CMAKE_CURRENT_BINARY_DIR}/, +# keeping generated .o/.ali files in the build tree (gnatmake -D) rather +# than polluting the source directory. The custom target is named +# ada_, not : target names must be unique across the whole +# combined root project (every PRK language directory has its own p2p/ +# stencil/etc.), same reasoning as every other language's CMakeLists.txt +# in this repo -- but since this is a custom target/command, not a real +# add_executable(), there's no OUTPUT_NAME to reconcile: the produced +# binary is simply named directly. +function(prk_ada_executable name) + set(_src "${CMAKE_CURRENT_LIST_DIR}/${name}.adb") + set(_out "${CMAKE_CURRENT_BINARY_DIR}/${name}") + add_custom_command( + OUTPUT "${_out}" + COMMAND "${GNATMAKE_EXECUTABLE}" -D "${CMAKE_CURRENT_BINARY_DIR}" -p + "${_src}" -o "${_out}" + DEPENDS "${_src}" + WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" + COMMENT "Building Ada ${name} with gnatmake") + add_custom_target(ada_${name} ALL DEPENDS "${_out}") +endfunction() + +function(prk_ada_test test_name) + set(multiValueArgs COMMAND) + cmake_parse_arguments(PAT "" "" "${multiValueArgs}" ${ARGN}) + add_test(NAME ${test_name} COMMAND ${PAT_COMMAND}) + set_tests_properties(${test_name} PROPERTIES TIMEOUT 30) +endfunction() + +prk_ada_executable(nstream) +prk_ada_executable(p2p) +prk_ada_executable(stencil) +prk_ada_executable(transpose) + +# --- CTest: same invocations used to verify each kernel by hand --------- + +prk_ada_test(ada_nstream COMMAND "${CMAKE_CURRENT_BINARY_DIR}/nstream" 10 16000000) +prk_ada_test(ada_p2p COMMAND "${CMAKE_CURRENT_BINARY_DIR}/p2p" 10 1000 1000) +prk_ada_test(ada_transpose COMMAND "${CMAKE_CURRENT_BINARY_DIR}/transpose" 10 1000) +prk_ada_test(ada_stencil_star COMMAND "${CMAKE_CURRENT_BINARY_DIR}/stencil" 10 1000) +prk_ada_test(ada_stencil_grid COMMAND "${CMAKE_CURRENT_BINARY_DIR}/stencil" 10 1000 grid) +prk_ada_test(ada_stencil_star_r4 COMMAND "${CMAKE_CURRENT_BINARY_DIR}/stencil" 10 1000 star 4) diff --git a/ADA/Makefile b/ADA/Makefile index b18bb0a12..d85edc053 100644 --- a/ADA/Makefile +++ b/ADA/Makefile @@ -1,12 +1,12 @@ AC = gnat ACFLAGS = make -all: hello args nstream nstream_array transpose +all: nstream p2p stencil transpose %: %.adb $(AC) $(ACFLAGS) $< -o $@ clean: -rm -f *.i *.o *.ali - -rm -f hello args nstream nstream_array transpose + -rm -f nstream p2p stencil transpose diff --git a/ADA/args.adb b/ADA/args.adb deleted file mode 100644 index 492f0b65e..000000000 --- a/ADA/args.adb +++ /dev/null @@ -1,42 +0,0 @@ -with - Ada.Text_IO, - Ada.Integer_Text_IO, - Ada.Strings, - Ada.Strings.Bounded, - Ada.Command_line; - -use - Ada.Text_IO, - Ada.Integer_Text_IO, - Ada.Strings, - Ada.Strings.Bounded, - Ada.Command_line; - -procedure Args is - - package BS is new Ada.Strings.Bounded.Generic_Bounded_Length (Max => 15); - use BS; - - S : Bounded_String; - - A : Integer := 1; - I : Integer := 10; - N : Integer := 1_000_000; - -begin - - Put_Line("Args, World!"); - - --Put("Argument_Count="); - --Put(Item => Argument_Count, Width => 1); - Put_Line("Argument_Count=" & Argument_Count'Image); - - if Argument_Count > 0 then - Put_Line("Arg1=" & Argument(1)); - end if; - if Argument_Count > 1 then - Put_Line("Arg2=" & Argument(2)); - end if; - -end Args; - diff --git a/ADA/hello.adb b/ADA/hello.adb deleted file mode 100644 index 9356c17b8..000000000 --- a/ADA/hello.adb +++ /dev/null @@ -1,6 +0,0 @@ -with Text_IO; use Text_IO; -procedure Hello is -begin - Put_Line("Hello, World!"); -end Hello; - diff --git a/ADA/nstream.adb b/ADA/nstream.adb index 8096d8917..21051421c 100644 --- a/ADA/nstream.adb +++ b/ADA/nstream.adb @@ -12,8 +12,9 @@ procedure Nstream is Ada.Real_Time, Ada.Command_line; + -- GNAT Integer = int32_t, Long_Integer = int64_t Iterations : Integer := 10; - Length : Integer := 1_000_000; + Length : Long_Integer := 1_000_000; Scalar : constant := 3.0; @@ -26,23 +27,28 @@ begin Iterations := Integer'Value(Argument(1)); end if; if Argument_Count > 1 then - Length := Integer'Value(Argument(2)); + Length := Long_Integer'Value(Argument(2)); end if; - if Iterations < 2 then + if Iterations < 1 then Put_Line("Iteration count must be greater than " & Integer'Image(Iterations) ); end if; Put_Line("Number of iterations =" & Integer'Image(Iterations) ); - Put_Line("Vector length =" & Integer'Image(Length) ); + Put_Line("Vector length =" & Long_Integer'Image(Length) ); declare - type Vector is array(1..Length) of Float; - - I : Integer := 0; - A : Vector; - B : Vector; - C : Vector; + -- Heap-allocated (via an access type), not a plain local array: + -- a stack-declared array of this size does not reliably fit on + -- the default 8MB stack once Length grows past a few hundred + -- thousand elements. Double precision throughout, matching the + -- "double precision" convention of the other language ports. + type Long_Float_Array is array(Long_Integer range <>) of Long_Float + with Default_Component_Value => 0.0; + type Long_Float_Array_Access is access Long_Float_Array; + A : Long_Float_Array_Access := new Long_Float_Array(1..Length); + B : Long_Float_Array_Access := new Long_Float_Array(1..Length); + C : Long_Float_Array_Access := new Long_Float_Array(1..Length); T0, T1 : Time; DT : Time_Span; @@ -50,39 +56,38 @@ begin US : constant Time_Span := Microseconds(US => Iterations); K : Integer := 0; - AR : Float := 0.0; - BR : Float := 2.0; - CR : Float := 2.0; - Asum : Float := 0.0; + AR : Long_Float := 0.0; + BR : Long_Float := 2.0; + CR : Long_Float := 2.0; + Asum : Long_Float := 0.0; - AvgTime : Duration; - Bytes : Integer := Float'Size / 8; + Bytes : Long_Integer; Nstream_us : Integer; - Nstream_time : Float; - Bandwidth : Float; + Nstream_time : Long_Float; + Bandwidth : Long_Float; begin -- initialization for I in 1..Length Loop - A(I) := Float(0); - B(I) := Float(2); - C(I) := Float(2); + A(I) := 0.0; + B(I) := 2.0; + C(I) := 2.0; end Loop; - + -- run the experiment for K in 0..Iterations Loop - + if K = 1 then T0 := Clock; end if; - + for I in 1..Length Loop A(I) := A(I) + B(I) + Scalar * C(I); end Loop; - + end Loop; T1 := Clock; DT := T1 - T0; @@ -98,22 +103,18 @@ begin end Loop; if Asum /= 0.0 then - Put_Line("Asum=" & Float'Image(Asum) ); + Put_Line("Asum=" & Long_Float'Image(Asum) ); else Put_Line("Solution validates"); - Bytes := Bytes * Length * 4; + -- 4 array touches (read A,B,C + write A) * 8 bytes each + Bytes := Length * Long_Float'Size / 2; Nstream_us := DT / US; -- this is per iteration now, thanks to US - Nstream_time := Float(Nstream_us) / Float(1000000); - Put_Line("Avg time (s): " & Float'Image(Nstream_time)); - Bandwidth := Float(Bytes) / Float(Nstream_us); - Put_Line("Rate (MB/s): " & Float'Image(Bandwidth)); - -- archived for posterity - --Put_Line("Bytes=" & Integer'Image(Bytes) ); - AvgTime := To_Duration(DT); - Put_Line("Total Time: " & Duration'Image(AvgTime) & " seconds"); + Nstream_time := Long_Float(Nstream_us) / Long_Float(1_000_000); + Put_Line("Avg time (s): " & Long_Float'Image(Nstream_time)); + Bandwidth := Long_Float(Bytes) / Long_Float(Nstream_us); + Put_Line("Rate (MB/s): " & Long_Float'Image(Bandwidth)); end if; end; end Nstream; - diff --git a/ADA/nstream_array.adb b/ADA/nstream_array.adb deleted file mode 100644 index 52f5141a8..000000000 --- a/ADA/nstream_array.adb +++ /dev/null @@ -1,122 +0,0 @@ -with - Ada.Text_IO, - Ada.Integer_Text_IO, - Ada.Real_Time, - Ada.Command_line; - -procedure nstream_array is - - use - Ada.Text_IO, - Ada.Integer_Text_IO, - Ada.Real_Time, - Ada.Command_line; - - -- GNAT Integer = int32_t, Long_Integer = int64_t - Iterations : Integer := 10; - Length : Long_Integer := 1_000_000; - - Scalar : constant := 3.0; - -begin - - Put_Line("Parallel Research Kernels"); - Put_Line("Ada Serial STREAM triad: A = B + scalar * C"); - - if Argument_Count > 0 then - Iterations := Integer'Value(Argument(1)); - end if; - if Argument_Count > 1 then - Length := Long_Integer'Value(Argument(2)); - end if; - - if Iterations < 2 then - Put_Line("Iteration count must be greater than " & Integer'Image(Iterations) ); - end if; - - Put_Line("Number of iterations =" & Integer'Image(Iterations) ); - Put_Line("Vector length =" & Long_Integer'Image(Length) ); - - declare - type Long_Float_Array is array(Long_Integer Range <>) of Long_Float with Default_Component_Value => 0.0; - - I : Integer := 0; - A : access Long_Float_Array := new Long_Float_Array(1..Length); - B : access Long_Float_Array := new Long_Float_Array(1..Length); - C : access Long_Float_Array := new Long_Float_Array(1..Length); - - T0, T1 : Time; - DT : Time_Span; - - US : constant Time_Span := Microseconds(US => Iterations); - - K : Integer := 0; - AR : Long_Float := 0.0; - BR : Long_Float := 2.0; - CR : Long_Float := 2.0; - Asum : Long_Float := 0.0; - - AvgTime : Duration; - Bytes : Long_Integer; - Nstream_us : Integer; - Nstream_time : Long_Float; - Bandwidth : Long_Float; - - begin - --- initialization - - for I in 1..Length Loop - A(I) := Long_Float(0); - B(I) := Long_Float(2); - C(I) := Long_Float(2); - end Loop; - --- run the experiment - - for K in 0..Iterations Loop - - if K = 1 then - T0 := Clock; - end if; - - for I in 1..Length Loop - A(I) := A(I) + B(I) + Scalar * C(I); - end Loop; - - end Loop; - T1 := Clock; - DT := T1 - T0; - --- validation - - for K in 0..Iterations Loop - AR := AR + BR + Scalar * CR; - end Loop; - - for I in 1..Length Loop - Asum := Asum + ABS ( A(I) - AR ); - end Loop; - - if Asum /= 0.0 then - Put_Line("Asum=" & Long_Float'Image(Asum) ); - else - Put_Line("Solution validates"); - Bytes := Length * Long_Float'Size / 2; - Nstream_us := DT / US; -- this is per iteration now, thanks to US - Nstream_time := Long_Float(Nstream_us) / Long_Float(1000000); - Put_Line("Avg time (s): " & Long_Float'Image(Nstream_time)); - Bandwidth := Long_Float(Bytes) / Long_Float(Nstream_us); - Put_Line("Rate (MB/s): " & Long_Float'Image(Bandwidth)); - -- archived for posterity - --Put_Line("Bytes=" & Long_Integer'Image(Bytes) ); - --AvgTime := To_Duration(DT); - --Put_Line("Total Time: " & Duration'Image(AvgTime) & " seconds"); - --Put_Line("Integer'Last=" & Integer'Image(Integer'Last) ); - --Put_Line("Long_Integer'Last=" & Long_Integer'Image(Long_Integer'Last) ); - end if; - - end; - -end nstream_array; - diff --git a/ADA/p2p.adb b/ADA/p2p.adb new file mode 100644 index 000000000..48a504b35 --- /dev/null +++ b/ADA/p2p.adb @@ -0,0 +1,125 @@ +with + Ada.Text_IO, + Ada.Integer_Text_IO, + Ada.Real_Time, + Ada.Command_line; + +procedure p2p is + + use + Ada.Text_IO, + Ada.Integer_Text_IO, + Ada.Real_Time, + Ada.Command_line; + + -- GNAT Integer = int32_t, Long_Integer = int64_t + Iterations : Natural := 10; + M : Natural := 1_000; + N : Natural := 1_000; + +begin + + Put_Line("Parallel Research Kernels"); + Put_Line("Ada Serial pipeline execution on 2D grid"); + + if Argument_Count > 0 then + Iterations := Integer'Value(Argument(1)); + end if; + if Argument_Count > 1 then + M := Integer'Value(Argument(2)); + end if; + if Argument_Count > 2 then + N := Integer'Value(Argument(3)); + end if; + + if Iterations < 1 then + Put_Line("Iteration count must be greater than " & Integer'Image(Iterations) ); + end if; + + Put_Line("Number of iterations =" & Integer'Image(Iterations) ); + Put_Line("Grid sizes =" & Integer'Image(M) & " *" & Integer'Image(N) ); + + declare + + -- Double precision throughout, matching the "double precision" + -- convention of the other language ports of this kernel: + -- Ada.Numerics.Real_Arrays.Real_Matrix is single-precision Float, + -- which is not accurate enough for the corner-value recurrence + -- below once M/N (and thus the number of accumulated terms) grow. + -- Heap-allocated (like nstream_array.adb), not a plain local + -- array: an M-by-N Long_Float grid at realistic benchmark sizes + -- (e.g. 1000x1000 = 8MB) does not reliably fit on the default + -- 8MB stack. + type Grid_Matrix is array(Integer range <>, Integer range <>) of Long_Float; + type Grid_Matrix_Access is access Grid_Matrix; + Grid : Grid_Matrix_Access := new Grid_Matrix(1..M,1..N); + + T0, T1 : Time; + DT : Time_Span; + + US : constant Time_Span := Microseconds(US => Iterations); + + K : Integer := 0; + + CornerVal : Long_Float; + Epsilon : constant Long_Float := 1.0E-8; + + AvgTime : Long_Float; + Nstream_us : Integer; + + begin + +-- initialization: left column and top row carry a boundary ramp, everything +-- else starts at zero, same as the other language ports' p2p kernel. + + for J in 1..N Loop + Grid(1,J) := Long_Float(J-1); + end Loop; + for I in 1..M Loop + Grid(I,1) := Long_Float(I-1); + end Loop; + for I in 2..M Loop + for J in 2..N Loop + Grid(I,J) := 0.0; + end Loop; + end Loop; + +-- run the experiment + + for K in 0..Iterations Loop + + if K = 1 then + T0 := Clock; + end if; + + for I in 2..M Loop + for J in 2..N Loop + Grid(I,J) := Grid(I-1,J) + Grid(I,J-1) - Grid(I-1,J-1); + end Loop; + end Loop; + + -- copy top right corner value to bottom left corner to create + -- a dependency between successive iterations + Grid(1,1) := -Grid(M,N); + + end Loop; + T1 := Clock; + DT := T1 - T0; + +-- validation + + CornerVal := Long_Float((Iterations+1) * (M+N-2)); + + if abs ( Grid(M,N) - CornerVal ) / CornerVal > Epsilon then + Put_Line("ERROR: checksum " & Long_Float'Image(Grid(M,N)) & " does not match verification value " & Long_Float'Image(CornerVal) ); + else + Put_Line("Solution validates"); + Nstream_us := DT / US; -- average time per iteration, in microseconds, thanks to US + AvgTime := Long_Float(Nstream_us) / Long_Float(1_000_000); + Put_Line("Avg time (s): " & Long_Float'Image(AvgTime) ); + Put_Line("Rate (MFlops/s): " & Long_Float'Image(1.0E-6 * 2.0 * Long_Float((M-1)*(N-1)) / AvgTime) ); + end if; + + end; + +end p2p; diff --git a/ADA/stencil.adb b/ADA/stencil.adb new file mode 100644 index 000000000..45130d513 --- /dev/null +++ b/ADA/stencil.adb @@ -0,0 +1,192 @@ +with + Ada.Text_IO, + Ada.Integer_Text_IO, + Ada.Real_Time, + Ada.Command_line; + +procedure stencil is + + use + Ada.Text_IO, + Ada.Integer_Text_IO, + Ada.Real_Time, + Ada.Command_line; + + -- GNAT Integer = int32_t, Long_Integer = int64_t + Iterations : Natural := 10; + N : Natural := 1_000; + Radius : Natural := 2; + Is_Star : Boolean := True; + +begin + + Put_Line("Parallel Research Kernels"); + Put_Line("Ada Serial stencil execution on 2D grid"); + + if Argument_Count > 0 then + Iterations := Integer'Value(Argument(1)); + end if; + if Argument_Count > 1 then + N := Integer'Value(Argument(2)); + end if; + if Argument_Count > 2 then + Is_Star := Argument(3) = "star"; + end if; + if Argument_Count > 3 then + Radius := Integer'Value(Argument(4)); + end if; + + if Iterations < 1 then + Put_Line("Iteration count must be greater than " & Integer'Image(Iterations) ); + end if; + if 2*Radius + 1 > N then + Put_Line("ERROR: Stencil radius exceeds grid size"); + return; + end if; + + Put_Line("Number of iterations =" & Integer'Image(Iterations) ); + Put_Line("Grid size =" & Integer'Image(N) ); + if Is_Star then + Put_Line("Type of stencil = star"); + else + Put_Line("Type of stencil = grid"); + end if; + Put_Line("Radius of stencil =" & Integer'Image(Radius) ); + + declare + + -- Double precision throughout: Ada.Numerics.Real_Arrays.Real_Matrix + -- is single-precision Float, which does not have enough mantissa + -- bits to hold the L1-norm accumulation below to 1.0E-8 once the + -- stencil has more than a handful of points (grid pattern, or a + -- larger radius). + -- Weight is small (a (2*Radius+1)^2 stencil) and stays on the + -- stack; Input/Output are heap-allocated (like nstream_array.adb) + -- since an N-by-N Long_Float grid at realistic benchmark sizes + -- (e.g. 1000x1000 = 8MB) does not reliably fit on the default + -- 8MB stack. + type Real_Matrix is array(Integer range <>, Integer range <>) of Long_Float; + type Real_Matrix_Access is access Real_Matrix; + Weight : Real_Matrix(-Radius..Radius,-Radius..Radius) := (others => (others => 0.0)); + Input : Real_Matrix_Access := new Real_Matrix(0..N-1,0..N-1); + Output : Real_Matrix_Access := new Real_Matrix(0..N-1,0..N-1); + + T0, T1 : Time; + DT : Time_Span; + + US : constant Time_Span := Microseconds(US => Iterations); + + K : Integer := 0; + + StencilSize : Integer; + Norm : Long_Float := 0.0; + ActivePoints : constant Long_Float := Long_Float((N - 2*Radius) * (N - 2*Radius)); + ReferenceNorm : constant Long_Float := Long_Float(2 * (Iterations+1)); + Epsilon : constant Long_Float := 1.0E-8; + + AvgTime : Long_Float; + Nstream_us : Integer; + + begin + +-- initialize the input array: in(i,j) = i+j, same convention as the other +-- language ports of this kernel + + for I in 0..N-1 Loop + for J in 0..N-1 Loop + Input(I,J) := Long_Float(I+J); + Output(I,J) := 0.0; + end Loop; + end Loop; + +-- fill the weight stencil: a discrete approximation to the derivative +-- operator, same construction as the other language ports + + if Is_Star then + StencilSize := 4*Radius + 1; + for I in 1..Radius Loop + Weight(0,I) := 1.0 / Long_Float(2*I*Radius); + Weight(I,0) := 1.0 / Long_Float(2*I*Radius); + Weight(0,-I) := -1.0 / Long_Float(2*I*Radius); + Weight(-I,0) := -1.0 / Long_Float(2*I*Radius); + end Loop; + else + StencilSize := (2*Radius+1) * (2*Radius+1); + for J in 1..Radius Loop + for I in -J+1..J-1 Loop + Weight(I,J) := 1.0 / Long_Float(4*J*(2*J-1)*Radius); + Weight(I,-J) := -1.0 / Long_Float(4*J*(2*J-1)*Radius); + Weight(J,I) := 1.0 / Long_Float(4*J*(2*J-1)*Radius); + Weight(-J,I) := -1.0 / Long_Float(4*J*(2*J-1)*Radius); + end Loop; + Weight(J,J) := 1.0 / Long_Float(4*J*Radius); + Weight(-J,-J) := -1.0 / Long_Float(4*J*Radius); + end Loop; + end if; + +-- run the experiment + + for K in 0..Iterations Loop + + if K = 1 then + T0 := Clock; + end if; + + -- apply the stencil operator + for I in Radius..N-1-Radius Loop + for J in Radius..N-1-Radius Loop + if Is_Star then + for JJ in -Radius..Radius Loop + Output(I,J) := Output(I,J) + Weight(0,JJ) * Input(I,J+JJ); + end Loop; + for II in -Radius..-1 Loop + Output(I,J) := Output(I,J) + Weight(II,0) * Input(I+II,J); + end Loop; + for II in 1..Radius Loop + Output(I,J) := Output(I,J) + Weight(II,0) * Input(I+II,J); + end Loop; + else + for JJ in -Radius..Radius Loop + for II in -Radius..Radius Loop + Output(I,J) := Output(I,J) + Weight(II,JJ) * Input(I+II,J+JJ); + end Loop; + end Loop; + end if; + end Loop; + end Loop; + + -- add a constant to the input array to force a refresh of the + -- neighbor data on the next iteration + for I in 0..N-1 Loop + for J in 0..N-1 Loop + Input(I,J) := Input(I,J) + 1.0; + end Loop; + end Loop; + + end Loop; + T1 := Clock; + DT := T1 - T0; + +-- validation: L1 norm of the interior of the output array should match a +-- known closed-form reference value + + for I in Radius..N-1-Radius Loop + for J in Radius..N-1-Radius Loop + Norm := Norm + abs ( Output(I,J) ); + end Loop; + end Loop; + Norm := Norm / ActivePoints; + + if abs ( Norm - ReferenceNorm ) > Epsilon then + Put_Line("ERROR: L1 norm = " & Long_Float'Image(Norm) & " Reference L1 norm = " & Long_Float'Image(ReferenceNorm) ); + else + Put_Line("Solution validates"); + Nstream_us := DT / US; -- average time per iteration, in microseconds, thanks to US + AvgTime := Long_Float(Nstream_us) / Long_Float(1_000_000); + Put_Line("Avg time (s): " & Long_Float'Image(AvgTime) ); + Put_Line("Rate (MFlops/s): " & Long_Float'Image(1.0E-6 * Long_Float(2*StencilSize+1) * ActivePoints / AvgTime) ); + end if; + + end; + +end stencil; diff --git a/ADA/transpose.adb b/ADA/transpose.adb index e0cd3093c..07515f1ce 100644 --- a/ADA/transpose.adb +++ b/ADA/transpose.adb @@ -1,5 +1,4 @@ with - Ada.Numerics.Real_Arrays, Ada.Text_IO, Ada.Integer_Text_IO, Ada.Real_Time, @@ -8,7 +7,6 @@ with procedure transpose is use - Ada.Numerics.Real_Arrays, Ada.Text_IO, Ada.Integer_Text_IO, Ada.Real_Time, @@ -21,7 +19,7 @@ procedure transpose is begin Put_Line("Parallel Research Kernels"); - Put_Line("Ada Serial Transpose B += A^T"); + Put_Line("Ada Serial Matrix transpose: B += A^T"); if Argument_Count > 0 then Iterations := Integer'Value(Argument(1)); @@ -30,7 +28,7 @@ begin Order := Integer'Value(Argument(2)); end if; - if Iterations < 2 then + if Iterations < 1 then Put_Line("Iteration count must be greater than " & Integer'Image(Iterations) ); end if; @@ -39,10 +37,17 @@ begin declare - I : Natural := 1; - J : Natural := 1; - A : Real_Matrix(1..Order,1..Order); - B : Real_Matrix(1..Order,1..Order); + -- Double precision throughout: Ada.Numerics.Real_Arrays.Real_Matrix + -- is single-precision Float, not accurate enough to validate at + -- 1.0E-8 once Order (and the running sums in the closed-form + -- check below) get large. Heap-allocated (like nstream_array.adb): + -- two Order-by-Order Long_Float matrices at realistic benchmark + -- sizes (e.g. 1000x1000 = 8MB each) do not reliably fit on the + -- default 8MB stack. + type Real_Matrix is array(Integer range <>, Integer range <>) of Long_Float; + type Real_Matrix_Access is access Real_Matrix; + A : Real_Matrix_Access := new Real_Matrix(1..Order,1..Order); + B : Real_Matrix_Access := new Real_Matrix(1..Order,1..Order); T0, T1 : Time; DT : Time_Span; @@ -50,68 +55,70 @@ begin US : constant Time_Span := Microseconds(US => Iterations); K : Integer := 0; - Asum : Float := 0.0; + AbsErr : Long_Float := 0.0; + Addit : Long_Float; - AvgTime : Duration; + AvgTime : Long_Float; Bytes : Long_Integer; Nstream_us : Integer; - Nstream_time : Long_Float; - Bandwidth : Long_Float; + Epsilon : constant Long_Float := 1.0E-8; begin --- initialization +-- initialization: A(i,j) = i*order + j (0-based), same convention as the +-- other language ports of this kernel for I in 1..Order Loop for J in 1..Order Loop - A(I,J) := Float(0); - B(I,J) := Float(2); + A(I,J) := Long_Float((I-1)*Order + (J-1)); + B(I,J) := 0.0; end Loop; end Loop; - + -- run the experiment for K in 0..Iterations Loop - + if K = 1 then T0 := Clock; end if; - + + -- transpose the matrix: B += A^T, then bump A so successive + -- iterations touch fresh data (matches every other language's + -- transpose kernel) for I in 1..Order Loop - A(I,J) := A(I,J) + B(I,J); + for J in 1..Order Loop + B(J,I) := B(J,I) + A(I,J); + A(I,J) := A(I,J) + 1.0; + end Loop; end Loop; - + end Loop; T1 := Clock; DT := T1 - T0; --- validation +-- validation: every B(j,i) should equal the closed-form sum of the +-- (Iterations+1) transposed-and-incremented values of A(i,j) - for I in 1..Order Loop - for J in 1..Order Loop - Asum := Asum + ABS ( A(I,J) ); + Addit := (Long_Float(Iterations+1) * Long_Float(Iterations)) / 2.0; + for J in 1..Order Loop + for I in 1..Order Loop + AbsErr := AbsErr + abs ( B(J,I) + - ( Long_Float((Order*(I-1) + (J-1))) * Long_Float(Iterations+1) + Addit ) ); end Loop; end Loop; - if Asum /= 0.0 then - Put_Line("Asum=" & Float'Image(Asum) ); + if AbsErr > Epsilon then + Put_Line("ERROR: Aggregate squared error " & Long_Float'Image(AbsErr) & " exceeds threshold " & Long_Float'Image(Epsilon) ); else Put_Line("Solution validates"); - Bytes := Long_Integer(Order) * Long_Integer(Order) * Long_Float'Size / 4; + Bytes := 2 * Long_Integer(Long_Float'Size / 8) * Long_Integer(Order) * Long_Integer(Order); Nstream_us := DT / US; -- this is per iteration now, thanks to US - Nstream_time := Long_Float(Nstream_us) / Long_Float(1000000); - Put_Line("Avg time (s): " & Long_Float'Image(Nstream_time)); - Bandwidth := Long_Float(Bytes) / Long_Float(Nstream_us); - Put_Line("Rate (MB/s): " & Long_Float'Image(Bandwidth)); - -- archived for posterity - --Put_Line("Bytes=" & Long_Integer'Image(Bytes) ); - --AvgTime := To_Duration(DT); - --Put_Line("Total Time: " & Duration'Image(AvgTime) & " seconds"); - --Put_Line("Integer'Last=" & Integer'Image(Integer'Last) ); - --Put_Line("Long_Integer'Last=" & Long_Integer'Image(Long_Integer'Last) ); + AvgTime := Long_Float(Nstream_us) / Long_Float(1_000_000); + Put_Line("Avg time (s): " & Long_Float'Image(AvgTime) ); + Put_Line("Rate (MB/s): " & Long_Float'Image(Long_Float(Bytes) / Long_Float(Nstream_us)) ); end if; end; end transpose; - diff --git a/CMakeLists.txt b/CMakeLists.txt index c9bb51031..46704dcab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,6 +24,7 @@ option(PRK_ENABLE_MPIOPENMP "Build MPIOPENMP/ kernels" ON) option(PRK_ENABLE_MPIRMA "Build MPIRMA/ kernels" ON) option(PRK_ENABLE_MPISHM "Build MPISHM/ kernels" ON) option(PRK_ENABLE_SHMEM "Build SHMEM/ kernels" ON) +option(PRK_ENABLE_ADA "Build ADA/ kernels" ON) if(PRK_ENABLE_FORTRAN) add_subdirectory(FORTRAN) @@ -65,6 +66,10 @@ if(PRK_ENABLE_SHMEM) add_subdirectory(SHMEM) endif() +if(PRK_ENABLE_ADA) + add_subdirectory(ADA) +endif() + # Future languages follow the same pattern: add an option # PRK_ENABLE_ (default ON) and add_subdirectory(), with a # /CMakeLists.txt structured like FORTRAN/CMakeLists.txt -- its own