EPCD11006 Numerical Algorithms for High Performance Computing

Released: Monday 24th March 2025

Due: Monday 31st March 2025 by 11:59am (GMT)

In lectures we saw that an important way in which linear algebra libraries such  as LAPACK and ScaLAPACK obtain good performance is by somehow dividing  matrices up into blocks when executing algorithms such as LU factorisation.  This coursework explores this topic further. You are asked to conduct a practical investigation that requires you to apply knowledge gained from lectures and  practicals supplemented with limited consultation of external resources such  as the official (Sca)LAPACK documentation.

Marks are available for the correctness and quality of your code and scripts, the explanations you give in response to the questions, and the quality of the plots you are asked to make.  Indications of marks achievable for each part below constitute a total for all aspects of your submission related to that part.

Please ensure that you provide the following in your submission (using sensible file names throughout):

• Your code, which should compile and run on Cirrus as outlined below

• Your written answers to the questions below

• The plots you are asked to generate

• The raw data you used to generate these plots

• Any job scripts you used to run experiments

Written answers and plots should be given in one document, while code, data and job scripts should be provided in a separate zip or tar archive.

In addition to course materials you may want to consult (and reference if they directly inform. your answer) the following external resources:

• The User Guide for Netlib’s reference LAPACK and ScaLAPACK implementations, on which MKL was based:

https://www. netlib. org/lapack/lug https://netlib. org/scalapack/slug

•  Source code documentation for Netlib LAPACK routines such as SGETRF, available from https://www. netlib. org/lapack/explore-html

• For additional background, Chapter 10 of

https://link. springer. com/chapter/10.1007/3-540-36574-5_10

(accessible through Springer Institutional login via the University of Edinburgh).

1. LAPACK

You have been provided with a simple  LU factorisation program in both C (lufact-lapack. c) and Fortran (lufact-lapack. f90), similar to the code you worked with in exercises in the course.  You can choose to work in either of these two languages.  The provided code generates a random matrix A and right-hand side b using the matgen function, then performs an LU factorisation using the LAPACK sgetrf function before solving the problem for b with the sgetrs function.   The  right-hand side b  is chosen in matgen such that the solution x should be a column vector of ones: x = (1, . . . , 1)T . The code uses single precision throughout, and you will find that several matrix and vector helper functions are included.

To compile the provided code on Cirrus, first load the Intel compilers and MKL modules:

module  load  intel-20 .4/compilers

module  load  intel-20 .4/cmkl

Then to compile and link against LAPACK and BLAS provided by MKL, include the -mkl compiler option, i.e. for the two languages respectively:

ifort  -mkl  -o  lufact-lapack  lufact-lapack. f90

icc  -mkl  -o  lufact-lapack  lufact-lapack . c

(a) One of the main ways  LAPACK achieves good serial performance on modern processors with a multilevel cache hierarchy (L1, L2, L3) is by

dividing the full matrix A into blocks, i.e. submatrices, when executing an algorithm such as LU factorisation.  How and why does the block-based implementation of an algorithm like LU factorisation in LAPACK typically

yield higher serial performance compared to its implementation in LINPACK?

Your answer should include which aspects of theoretical algorithmic performance matter for performance given real-life hardware bottlenecks, and why.          [6 marks]

(b)  Block size can not be specified explicitly when calling LAPACK routines   and instead is automatically decided internally based on heuristics encoded in the LAPACK routine ILAENV.

Modify the provided code to interrogate ILAENV and determine the block size chosen by SGETRF when it is called by your program. Summarise your findings regarding how the block size changes going from matrix A of size 1x1 up to size 5000x5000. Why do you think the block size changes the way it does?                     [6 marks]

(c)  Insert calls to omp_get_wtime() before and after the call to SGETRF and use these to determine the time taken for LU factorisation (you will need to compile with OpenMP support).  Run experiments on Cirrus compute nodes to determine how this time changes as the matrix size increases across the range described above. Remember to load the same modules as during compilation, follow standard practice for running serial jobs on Cirrus as described in the Cirrus user guide, and ensure you follow good benchmarking practice, including:

•  Request exclusive node access (#SBATCH  --exclusive) and #SBATCH  --cpus-per-task=1 in your job script(s)

•  Set export  OMP_NUM_THREADS=1 to isolate  results from any potential effect of multithreading within MKL

•  Use srun  --cpu-bind=cores to launch your executable

Plot the time taken for SGETRF against the linear dimension N of the matrix A. Can you discern any effect of the changes in block size?                           [6 marks]

ScaLAPACK

You are also provided with Fortran code (lufact-scalapack. f90) that uses ScaLAPACK to run SGETRF in parallel using MPI and a block-cyclic processor grid decomposition as shown in lectures.  There should be no need to make any changes to the code itself and you do not need to be familiar with Fortran to complete this part of the coursework. The code contains several new routines and variables to allow the ScaLAPACK library to work but is otherwise basically the same as for the serial version. Small changes have been made to make the array sizes dynamic and there are a few small changes to the matgen routine.

Before compiling, make sure the following modules are loaded:

module  load  mpt

module  load  intel-20 .4/compilers

module  load  intel-20 .4/cmkl

A Makefile is provided, so to compile simply type make.

The executable takes the following input arguments:

./lufact matsize nprocrow nproccol blocksize

where matsize sets the size of the global matrix A (i.e. gives A the dimensions matsize  × matsize),  and the processor grid is set to have nprocrow rows of blocks and nproccol columns of blocks.   Finally blocksize sets the size measured in number of elements of the side of the blocks along both row and column dimensions, i.e. the blocks are square with dimensions blocksize  × blocksize (apart from any blocks at the boundaries of the matrix that may be cut off if they do  not fit exactly).   You  should  take care to ensure that nprocrow × nproccol = nprocs (the total number of processors you run on). As before, the code is designed to give a solution of xi  = 1.  Rather than print out the entire x array, the code prints out the value of |x − 1|, where 1 is a vector with each element equal to 1.

(d)  Submit the code to Cirrus using a batch script that launches lufact on 16 processes and with a 4x4 processor grid.  Remember to load the same modules as during compilation, follow standard practice for MPI-parallel jobs on Cirrus as described in the Cirrus user guide, and ensure you follow good benchmarking practice as outlined above.  Try running the code with blocksize=6400 using different block sizes (eg 8, 16, 32, 64, 128, 256). What is the optimal block size for this particular problem? Can you give likely reasons for the trend in performance when using different block sizes?                                                   [7 marks]