PTfS22: Assignment 1: Loop kernels, π reloaded (due May 13)

Opened: Wednesday, 4 May 2022, 12:00 AM

Due: Friday, 13 May 2022, 12:00 PM

Loop kernel benchmarking is a variant of microbenchmarking - running simple stuff that lets you get insight into the inner workings of a machine. The first step in this process is to get it right and to present the data correctly. This is what this task is all about.
Write a benchmark program that measures the performance in MFlop/s of the following two loop kernels (loop over i is implied):

(a) (20 credits) vector triad: a[i] = b[i] + c[i] * d[i]
(b) (20 credits) vector update: a[i] = s * a[i]

a, b, c, d are double-precision arrays of length N. s is a double-precision scalar that you should define as a compile-time constant:
```
const double s = 1.00000000001;
```
Allocate memory for all arrays on the heap, i.e. , using malloc() in C, new in C++, or allocate() in Fortran90. Do not forget to initialize all data elements with valid floating-point (FP) numbers. Using calloc() is not sufficient - see this blog post for an explanation.

Run your code using the following vector lengths (in elements): N=int(1.5 ^k), k = 8 ... 44. This will give you a decent resolution and equidistant points on the x axis (see below).

Perform the measurement on one core of the Fritz cluster for the given loop lengths (do not forget to fix the clock frequency to 2 GHz). For reasons of accuracy, make sure that the runtime of each kernel with each vector/matrix size is larger than 0.1 seconds by repeating the computation kernel in an outer loop. Maybe it is a good idea to dynamically adjust the number of repetitions depending on the runtime of the kernel:
```
            NITER=1;
            do {
              // time measurement
              wct_start = getTimeStamp();

              // repeat measurement often enough
              for(k=0; k<NITER; ++k) {
                // This is the benchmark loop
                for(i=0; i<N; ++i) {
                 // put loop body here: a[i] = ...
                }
                // end of benchmark loop
                if(a[N/2]<0.) printf("%lf",a[N/2]); // prevent compiler from eliminating loop
              }
              wct_end = getTimeStamp();
              NITER = NITER*2;
            } while (wct_end-wct_start<0.1); // at least 100 ms

            NITER = NITER/2;

            printf("Total walltime: %f, NITER: %d\n",wct_end-wct_start,NITER);
            
```
Make sure that the operations in the kernel actually get executed - compilers are smart! (This is what the bogus if statement is for: The compiler must not be able to determine the result of the condition at compile time. Example: if(a[N/2]<0.) - if all arrays are initialized with positive numbers, this condition is never true.) Use the standard compiler options -O3 -xHost -fno-alias. A compilable skeleton code can be found in ~ptfs100h/GettingStarted/scan_c.c.

Use your favorite graphics program (e.g., gnuplot or xmgrace or LibreOffice Calc or ...) to generate plots of the performance in MFlop/s vs. N. Choose a logarithmic scale on the x axis (think about why this is advisable). Always let the y axis start at zero. If you don't, the graph may be misleading and you will collect some very bad karma.

(c) (40 credits) Now use your code to benchmark the "recursive vector update" in the same way as above: a[i] = s * a[i-1]

Here you will explore two compilers: Use the "old" Intel compiler (icc/icpc/ifort for C/C++/Fortran) first, and then, in a second experiment, use the new compiler (icx,/icpx/ifx). Both groups of compilers are available once you load the intel module. Compare the results. Discuss the fundamental difference and try to think of an optimization the "new" compiler can do which the "old" one cannot.
π reloaded. (You will only fully appreciate this exercise after the lecture on Monday, May 9.) In Assignment 0 you studied a program that calculated the value of π via numerical integration (midpoint rule). In double precision, one iteration took almost exactly 2 cycles, and this is entirely determined by the divide instructions (everything else happens so fast that it can be hidden behind them). Re-run your code with the following alternative compiler options and report the number of cycles per iteration (use the new icx/icpx/ifx compilers again):

(a) -O3 -no-vec (no SIMD, forces scalar code) (5 credits)
(b) -O3 -xHost -qopt-zmm-usage=high (allows for AVX-512 instructions) (5 credits)

Knowing that the compiler produces AVX2 instructions (vector length 256 bits) with -xHost alone and actually uses an instruction that does four double-precision divides, what is your conclusion from the observed scalar and AVX-512 performance? (FYI: With AVX-512, a SIMD divide instruction is used that does eight divides). Calculate the inverse throughput (in cycles per instruction) of the scalar, AVX2, and AVX-512 double-precision divide instructions on this machine! (5 credits) Discuss the result. What is your conclusion? (5 credits)