Note, these notes contain additional material om optimization and parallelization. Parts of this material will be discussed this week.

Why resampling methods ?

Statistical analysis

Our simulations can be treated as computer experiments. This is particularly the case for Monte Carlo methods
The results can be analysed with the same statistical tools as we would use analysing experimental data.
As in all experiments, we are looking for expectation values and an estimate of how accurate they are, i.e., possible sources for errors.

Statistical analysis

As in other experiments, many numerical experiments have two classes of errors:
1. Statistical errors
2. Systematical errors
Statistical errors can be estimated using standard tools from statistics
Systematical errors are method specific and must be treated differently from case to case.

And why do we use such methods?

As you will see below, due to correlations between various measurements, we need to evaluate the so-called covariance in order to establish a proper evaluation of the total variance and the thereby the standard deviation of a given expectation value.

The covariance however, leads to an evaluation of a double sum over the various stochastic variables. This becomes computationally too expensive to evaluate. Methods like the Bootstrap, the Jackknife and/or Blocking allow us to circumvent this problem.

Central limit theorem

Measurement $ i $

We assumed that each individual measurement $ x_{ij} $ is represented by stochastic variables which independent and identically distributed (iid). This defined the sample mean of of experiment $ i $ with $ n $ samples as

$$ \overline{x}_i=\frac{1}{n}\sum_{j} x_{ij}. $$

and the sample variance

$$ \sigma^2_i=\frac{1}{n}\sum_{j} \left(x_{ij}-\overline{x}_i\right)^2. $$

Further remarks

Note that we use $ n $ instead of $ n-1 $ in the definition of variance. The sample variance and the sample mean are not necessarily equal to the exact values we would get if we knew the corresponding probability distribution.

Running many measurements

Adding $ m $ measurements $ i $

With the assumption that the average measurements $ i $ are also defined as iid stochastic variables and have the same probability function $ p $, we defined the total average over $ m $ experiments as

$$ \overline{X}=\frac{1}{m}\sum_{i} \overline{x}_{i}. $$

and the total variance

$$ \sigma^2_{m}=\frac{1}{m}\sum_{i} \left( \overline{x}_{i}-\overline{X}\right)^2. $$

These are the quantities we used in showing that if the individual mean values are iid stochastic variables, then in the limit $ m\rightarrow \infty $, the distribution for $ \overline{X} $ is given by a Gaussian distribution with variance $ \sigma^2_m $.

Adding more definitions

and introducing the centered value $ \tilde{x}_{ij}=x_{ij}-\overline{X} $, we can rewrite $ \sigma_m^2 $ as

Further rewriting

We can rewrite the latter in terms of a sum over diagonal elements only and another sum which contains the non-diagonal elements

The first term on the last rhs is nothing but the total sample variance $ \sigma^2 $ divided by $ m $. The second term represents the covariance.

The covariance term

The first term is what we have used till now in order to estimate the standard deviation. However, the second term which gives us a measure of the correlations between different stochastic events, can result in contributions which give rise to a larger standard deviation and variance $ \sigma_m^2 $. Note also the evaluation of the second term leads to a double sum over all events. If we run a VMC calculation with say $ 10^9 $ Monte carlo samples, the latter term would lead to $ 10^{18} $ function evaluations. We don't want to, by obvious reasons, to venture into that many evaluations.

Note also that if our stochastic events are iid then the covariance terms is zero.

Rewriting the covariance term

Introducing the correlation function

We introduce then a correlation function $ \kappa_d=f_d/\sigma^2 $. Note that $ \kappa_0 =1 $. We rewrite the variance $ \sigma_m^2 $ as

The code here shows the evolution of $ \kappa_d $ as a function of $ d $ for a series of random numbers. We see that the function $ \kappa_d $ approaches $ 0 $ as $ d\rightarrow \infty $.

In this case, our data are given by random numbers generated for the uniform distribution with $ x\in [0,1] $. Even with two random numbers being far away, we note that the correlation function is not zero.

Computing the correlation function

Resampling methods: Blocking

The blocking method was made popular by Flyvbjerg and Pedersen (1989) and has become one of the standard ways to estimate the variance $ \mathrm{var}(\widehat{\theta}) $ for exactly one estimator $ \widehat{\theta} $, namely $ \widehat{\theta} = \overline{X} $, the mean value.

Assume $ n = 2^d $ for some integer $ d>1 $ and $ X_1,X_2,\cdots, X_n $ is a stationary time series to begin with. Moreover, assume that the series is asymptotically uncorrelated. We switch to vector notation by arranging $ X_1,X_2,\cdots,X_n $ in an $ n $-tuple. Define:

Why blocking?

The strength of the blocking method is when the number of observations, $ n $ is large. For large $ n $, the complexity of dependent bootstrapping scales poorly, but the blocking method does not, moreover, it becomes more accurate the larger $ n $ is.

Blocking Transformations

We now define the blocking transformations. The idea is to take the mean of subsequent pair of elements from $ \boldsymbol{X} $ and form a new vector $ \boldsymbol{X}_1 $. Continuing in the same way by taking the mean of subsequent pairs of elements of $ \boldsymbol{X}_1 $ we obtain $ \boldsymbol{X}_2 $, and so on. Define $ \boldsymbol{X}_i $ recursively by:

Blocking transformations

The quantity $ \boldsymbol{X}_k $ is subject to $ k $ blocking transformations. We now have $ d $ vectors $ \boldsymbol{X}_0, \boldsymbol{X}_1,\cdots,\boldsymbol{X}_{d-1} $ containing the subsequent averages of observations. It turns out that if the components of $ \boldsymbol{X} $ is a stationary time series, then the components of $ \boldsymbol{X}_i $ is a stationary time series for all $ 0 \leq i \leq d-1 $

We can then compute the autocovariance (or just covariance), the variance, sample mean, and number of observations for each $ i $. Let $ \gamma_i, \sigma_i^2, \overline{X}_i $ denote the covariance, variance and average of the elements of $ \boldsymbol{X}_i $ and let $ n_i $ be the number of elements of $ \boldsymbol{X}_i $. It follows by induction that $ n_i = n/2^i $.

Blocking Transformations

Using the definition of the blocking transformation and the distributive property of the covariance, it is clear that since $ h =|i-j| $ we can define

The quantity $ \hat{X} $ is asymptotically uncorrelated by assumption, $ \hat{X}_k $ is also asymptotic uncorrelated. Let's turn our attention to the variance of the sample mean $ \mathrm{var}(\overline{X}) $.

Blocking Transformations, getting there

We can show that $ \mathrm{var}(\overline{X}_i) = \mathrm{var}(\overline{X}_j) $ for all $ 0 \leq i \leq d-1 $ and $ 0 \leq j \leq d-1 $.

Blocking Transformations, final expressions

By repeated use of this equation we get $ \mathrm{var}(\overline{X}_i) = \mathrm{var}(\overline{X}_0) = \mathrm{var}(\overline{X}) $ for all $ 0 \leq i \leq d-1 $. This has the consequence that

More on optimization

We have additional compiler options for optimization. These may include procedure inlining where performance may be improved, moving constants inside loops outside the loop, identify potential parallelism, include automatic vectorization or replace a division with a reciprocal and a multiplication if this speeds up the code.

c++  -O3 -c  mycode.cpp
c++  -O3 -o  mycode.exe  mycode.o

This (other options are -O2 or -Ofast) is the recommended option.

Optimization and profiling

It is also useful to profile your program under the development stage. You would then compile with

c++  -pg -O3 -c  mycode.cpp
c++  -pg -O3 -o  mycode.exe  mycode.o

After you have run the code you can obtain the profiling information via

gprof mycode.exe >  ProfileOutput

When you have profiled properly your code, you must take out this option as it slows down performance. For memory tests use valgrind. An excellent environment for all these aspects, and much more, is Qt creator.

Optimization and debugging

Adding debugging options is a very useful alternative under the development stage of a program. You would then compile with

c++  -g -O0 -c  mycode.cpp
c++  -g -O0 -o  mycode.exe  mycode.o

This option generates debugging information allowing you to trace for example if an array is properly allocated. Some compilers work best with the no optimization option -O0.

Other optimization flags

Depending on the compiler, one can add flags which generate code that catches integer overflow errors. The flag -ftrapv does this for the CLANG compiler on OS X operating systems.

Other hints

In general, irrespective of compiler options, it is useful to

avoid if tests or call to functions inside loops, if possible.
avoid multiplication with constants inside loops if possible

Here is an example of a part of a program where specific operations lead to a slower code

k = n-1;
for (i = 0; i < n; i++){
    a[i] = b[i] +c*d;
    e = g[k];
}

A better code is

temp = c*d;
for (i = 0; i < n; i++){
    a[i] = b[i] + temp;
}
e = g[n-1];

Here we avoid a repeated multiplication inside a loop. Most compilers, depending on compiler flags, identify and optimize such bottlenecks on their own, without requiring any particular action by the programmer. However, it is always useful to single out and avoid code examples like the first one discussed here.

Vectorization and the basic idea behind parallel computing

Present CPUs are highly parallel processors with varying levels of parallelism. The typical situation can be described via the following three statements.

Pursuit of shorter computation time and larger simulation size gives rise to parallel computing.
Multiple processors are involved to solve a global problem.
The essence is to divide the entire computation evenly among collaborative processors. Divide and conquer.

Before we proceed with a more detailed discussion of topics like vectorization and parallelization, we need to remind ourselves about some basic features of different hardware models.

A rough classification of hardware models

Conventional single-processor computers are named SISD (single-instruction-single-data) machines.
SIMD (single-instruction-multiple-data) machines incorporate the idea of parallel processing, using a large number of processing units to execute the same instruction on different data.
Modern parallel computers are so-called MIMD (multiple-instruction-multiple-data) machines and can execute different instruction streams in parallel on different data.

Shared memory and distributed memory

One way of categorizing modern parallel computers is to look at the memory configuration.

In shared memory systems the CPUs share the same address space. Any CPU can access any data in the global memory.
In distributed memory systems each CPU has its own memory.

The CPUs are connected by some network and may exchange messages.

Different parallel programming paradigms

Task parallelism: the work of a global problem can be divided into a number of independent tasks, which rarely need to synchronize. Monte Carlo simulations represent a typical situation. Integration is another. However this paradigm is of limited use.
Data parallelism: use of multiple threads (e.g. one or more threads per processor) to dissect loops over arrays etc. Communication and synchronization between processors are often hidden, thus easy to program. However, the user surrenders much control to a specialized compiler. Examples of data parallelism are compiler-based parallelization and OpenMP directives.

Different parallel programming paradigms

Message passing: all involved processors have an independent memory address space. The user is responsible for partitioning the data/work of a global problem and distributing the subproblems to the processors. Collaboration between processors is achieved by explicit message passing, which is used for data transfer plus synchronization.
This paradigm is the most general one where the user has full control. Better parallel efficiency is usually achieved by explicit message passing. However, message-passing programming is more difficult.

What is vectorization?

Vectorization is a special case of Single Instructions Multiple Data (SIMD) to denote a single instruction stream capable of operating on multiple data elements in parallel. We can think of vectorization as the unrolling of loops accompanied with SIMD instructions.

Vectorization is the process of converting an algorithm that performs scalar operations (typically one operation at the time) to vector operations where a single operation can refer to many simultaneous operations. Consider the following example

If the code is not vectorized, the compiler will simply start with the first element and then perform subsequent additions operating on one address in memory at the time.

Number of elements that can acted upon

A SIMD instruction can operate on multiple data elements in one single instruction. It uses the so-called 128-bit SIMD floating-point register. In this sense, vectorization adds some form of parallelism since one instruction is applied to many parts of say a vector.

The number of elements which can be operated on in parallel range from four single-precision floating point data elements in so-called Streaming SIMD Extensions and two double-precision floating-point data elements in Streaming SIMD Extensions 2 to sixteen byte operations in a 128-bit register in Streaming SIMD Extensions 2. Thus, vector-length ranges from 2 to 16, depending on the instruction extensions used and on the data type.

Number of elements that can acted upon, examples

If the code is not vectorized and we have a 128-bit register to store a 32 bits floating point number, it means that we have $ 3\times 32 $ bits that are not used.

We have thus unused space in our SIMD registers. These registers could hold three additional integers.

Operation counts for scalar operation

Number of elements that can acted upon, examples

If we vectorize the code, we can perform, with a 128-bit register four simultaneous operations, that is we have

Number of operations when vectorized

A simple test case with and without vectorization

We implement these operations in a simple c++ program that computes at the end the norm of a vector.

Compiling with and without vectorization

The speedup depends on the size of the vectors. In the example here we have run with $ 10^7 $ elements. The example here was run on an IMac17.1 with OSX El Capitan (10.11.4) as operating system and an Intel i5 3.3 GHz CPU.

This particular C++ compiler speeds up the above loop operations with a factor of 1.5 Performing the same operations for $ 10^9 $ elements results in a smaller speedup since reading from main memory is required. The non-vectorized code is seemingly faster.

Compiling with and without vectorization using clang

Automatic vectorization and vectorization inhibitors, criteria

An important criteria is that the loop counter $ n $ is known at the entry of the loop.

The variable $ n $ does need to be known at compile time. However, this variable must stay the same for the entire duration of the loop. It implies that an exit statement inside the loop cannot be data dependent.

Automatic vectorization and vectorization inhibitors, exit criteria

An exit statement should in general be avoided. If the exit statement contains data-dependent conditions, the loop cannot be vectorized. The following is an example of a non-vectorizable loop

Avoid loop termination conditions and opt for a single entry loop variable $ n $. The lower and upper bounds have to be kept fixed within the loop.

Automatic vectorization and vectorization inhibitors, straight-line code

SIMD instructions perform the same type of operations multiple times. A switch statement leads thus to a non-vectorizable loop since different statemens cannot branch. The following code can however be vectorized since the if statement is implemented as a masked assignment.

These operations can be performed for all data elements but only those elements which the mask evaluates as true are stored. In general, one should avoid branches such as switch, go to, or return statements or if constructs that cannot be treated as masked assignments.

Automatic vectorization and vectorization inhibitors, nested loops

The exception is if an original outer loop is transformed into an inner loop as the result of compiler optimizations.

Automatic vectorization and vectorization inhibitors, function calls

Calls to programmer defined functions ruin vectorization. However, calls to intrinsic functions like $ \sin{x} $, $ \cos{x} $, $ \exp{x} $ etc are allowed since they are normally efficiently vectorized. The following example is fully vectorizable

Similarly, inline functions defined by the programmer, allow for vectorization since the function statements are glued into the actual place where the function is called.

Automatic vectorization and vectorization inhibitors, data dependencies

One has to keep in mind that vectorization changes the order of operations inside a loop. A so-called read-after-write statement with an explicit flow dependency cannot be vectorized. The following code

is an example of flow dependency and results in wrong numerical results if vectorized. For a scalar operation, the value $ a[i-1] $ computed during the iteration is loaded into the right-hand side and the results are fine. In vector mode however, with a vector length of four, the values $ a[0] $, $ a[1] $, $ a[2] $ and $ a[3] $ from the previous loop will be loaded into the right-hand side and produce wrong results. That is, we have

and if the two first iterations are executed at the same by the SIMD instruction, the value of say $ a[1] $ could be used by the second iteration before it has been calculated by the first iteration, leading thereby to wrong results.

Automatic vectorization and vectorization inhibitors, more data dependencies

On the other hand, a so-called write-after-read statement can be vectorized. The following code

is an example of flow dependency that can be vectorized since no iteration with a higher value of $ i $ can complete before an iteration with a lower value of $ i $. However, such code leads to problems with parallelization.

Automatic vectorization and vectorization inhibitors, memory stride

For C++ programmers it is also worth keeping in mind that an array notation is preferred to the more compact use of pointers to access array elements. The compiler can often not tell if it is safe to vectorize the code.

When dealing with arrays, you should also avoid memory stride, since this slows down considerably vectorization. When you access array element, write for example the inner loop to vectorize using unit stride, that is, access successively the next array element in memory, as shown here

Memory management

Loads and stores to memory can be as important as floating point operations when we measure performance.

Memory and communication

Many of these performance features are not captured in most programming languages.

Measuring performance

Problems with measuring time

Problems with cold start

What happens when the code is executed? The assumption is that the code is ready to execute. But

Problems with smart compilers

Problems with interference

Problems with measuring performance

Thomas algorithm for tridiagonal linear algebra equations

$$ \left( \begin{array}{ccccc} b_0 & c_0 & & & \\ a_0 & b_1 & c_1 & & \\ & & \ddots & & \\ & & a_{m-3} & b_{m-2} & c_{m-2} \\ & & & a_{m-2} & b_{m-1} \end{array} \right) \left( \begin{array}{c} x_0 \\ x_1 \\ \vdots \\ x_{m-2} \\ x_{m-1} \end{array} \right)=\left( \begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_{m-2} \\ f_{m-1} \\ \end{array} \right) $$

Thomas algorithm, forward substitution

The first step is to multiply the first row by $ a_0/b_0 $ and subtract it from the second row. This is known as the forward substitution step. We obtain then

$$ a_i = 0, $$ $$ b_i = b_i - \frac{a_{i-1}}{b_{i-1}}c_{i-1}, $$

and

$$ f_i = f_i - \frac{a_{i-1}}{b_{i-1}}f_{i-1}. $$

At this point the simplified equation, with only an upper triangular matrix takes the form

$$ \left( \begin{array}{ccccc} b_0 & c_0 & & & \\ & b_1 & c_1 & & \\ & & \ddots & & \\ & & & b_{m-2} & c_{m-2} \\ & & & & b_{m-1} \end{array} \right)\left( \begin{array}{c} x_0 \\ x_1 \\ \vdots \\ x_{m-2} \\ x_{m-1} \end{array} \right)=\left( \begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_{m-2} \\ f_{m-1} \\ \end{array} \right) $$

Thomas algorithm, backward substitution

The next step is the backward substitution step. The last row is multiplied by $ c_{N-3}/b_{N-2} $ and subtracted from the second to last row, thus eliminating $ c_{N-3} $ from the last row. The general backward substitution procedure is

$$ c_i = 0, $$

and

$$ f_{i-1} = f_{i-1} - \frac{c_{i-1}}{b_i}f_i $$

All that ramains to be computed is the solution, which is the very straight forward process of

$$ x_i = \frac{f_i}{b_i} $$

Thomas algorithm and counting of operations (floating point and memory)

We have in specific case the following operations with the floating operations

Memory Reads: $ 14(N-2) $;
Memory Writes: $ 4(N-2) $;
Subtractions: $ 3(N-2) $;
Multiplications: $ 3(N-2) $;
Divisions: $ 4(N-2) $.

// Forward substitution    
// Note that we can simplify by precalculating a[i-1]/b[i-1]
  for (int i=1; i < n; i++) {
     b[i] = b[i] - (a[i-1]*c[i-1])/b[i-1];
     f[i] = g[i] - (a[i-1]*f[i-1])/b[i-1];
  }
  x[n-1] = f[n-1] / b[n-1];
  // Backwards substitution                                                           
  for (int i = n-2; i >= 0; i--) {
     f[i] = f[i] - c[i]*f[i+1]/b[i+1];
     x[i] = f[i]/b[i];
  }

Example: Transpose of a matrix

Matrix-matrix multiplication

This the matrix-matrix multiplication code with plain c++ memory allocation. It computes at the end the Frobenius norm.

How do we define speedup? Simplest form

Speedup measures the ratio of performance between two objects
Versions of same code, with different number of processors
Serial and vector versions
Try different programing languages, c++ and Fortran
Two algorithms computing the same result

How do we define speedup? Correct baseline

The key is choosing the correct baseline for comparison

For our serial vs. vectorization examples, using compiler-provided vectorization, the baseline is simple; the same code, with vectorization turned off

For parallel applications, this is much harder:

Choice of algorithm, decomposition, performance of baseline case etc.

Parallel speedup

For parallel applications, speedup is typically defined as

Speedup $ =T_1/T_p $

Here $ T_1 $ is the time on one processor and $ T_p $ is the time using $ p $ processors.

Can the speedup become larger than $ p $? That means using $ p $ processors is more than $ p $ times faster than using one processor.

Speedup and memory

The speedup on $ p $ processors can be greater than $ p $ if memory usage is optimal! Consider the case of a memorybound computation with $ M $ words of memory

If $ M/p $ fits into cache while $ M $ does not, the time to access memory will be different in the two cases:
$ T_1 $ uses the main memory bandwidth
$ T_p $ uses the appropriate cache bandwidth

Upper bounds on speedup

Assume that almost all parts of a code are perfectly parallelizable (fraction $ f $). The remainder, fraction $ (1-f) $ cannot be parallelized at all.

That is, there is work that takes time $ W $ on one process; a fraction $ f $ of that work will take time $ Wf/p $ on $ p $ processors.

What is the maximum possible speedup as a function of $ f $?

Amdahl's law

On one processor we have

$$ T_1 = (1-f)W + fW = W $$

On $ p $ processors we have

$$ T_p = (1-f)W + \frac{fW}{p}, $$

resulting in a speedup of

$$ \frac{T_1}{T_p} = \frac{W}{(1-f)W+fW/p} $$

As $ p $ goes to infinity, $ fW/p $ goes to zero, and the maximum speedup is

$$ \frac{1}{1-f}, $$

meaning that if if $ f = 0.99 $ (all but $ 1\% $ parallelizable), the maximum speedup is $ 1/(1-.99)=100 $!

How much is parallelizable

If any non-parallel code slips into the application, the parallel performance is limited.

In many simulations, however, the fraction of non-parallelizable work is $ 10^{-6} $ or less due to large arrays or objects that are perfectly parallelizable.

Today's situation of parallel computing

Distributed memory is the dominant hardware configuration. There is a large diversity in these machines, from MPP (massively parallel processing) systems to clusters of off-the-shelf PCs, which are very cost-effective.
Message-passing is a mature programming paradigm and widely accepted. It often provides an efficient match to the hardware. It is primarily used for the distributed memory systems, but can also be used on shared memory systems.
Modern nodes have nowadays several cores, which makes it interesting to use both shared memory (the given node) and distributed memory (several nodes with communication). This leads often to codes which use both MPI and OpenMP.

Our lectures will focus on both MPI and OpenMP.

Overhead present in parallel computing

Uneven load balance: not all the processors can perform useful work at all time.
Overhead of synchronization
Overhead of communication
Extra computation due to parallelization

Due to the above overhead and that certain parts of a sequential algorithm cannot be parallelized we may not achieve an optimal parallelization.

Parallelizing a sequential algorithm

Identify the part(s) of a sequential algorithm that can be executed in parallel. This is the difficult part,
Distribute the global work and data among $ P $ processors.

Strategies

Develop codes locally, run with some few processes and test your codes. Do benchmarking, timing and so forth on local nodes, for example your laptop or PC.
When you are convinced that your codes run correctly, you can start your production runs on available supercomputers.

How do I run MPI on a PC/Laptop? MPI

To install MPI is rather easy on hardware running unix/linux as operating systems, follow simply the instructions from the OpenMPI website. See also subsequent slides. When you have made sure you have installed MPI on your PC/laptop,

Compile with mpicxx/mpic++ or mpif90

  # Compile and link
  mpic++ -O3 -o nameofprog.x nameofprog.cpp
  #  run code with for example 8 processes using mpirun/mpiexec
  mpiexec -n 8 ./nameofprog.x

Can I do it on my own PC/laptop? OpenMP installation

If you wish to install MPI and OpenMP on your laptop/PC, we recommend the following:

For OpenMP, the compile option -fopenmp is included automatically in recent versions of the C++ compiler and Fortran compilers. For users of different Linux distributions, simply use the available C++ or Fortran compilers and add the above compiler instructions, see also code examples below.
For OS X users however, install libomp

  brew install libomp

and compile and link as

c++ -o <name executable> <name program.cpp>  -lomp

Installing MPI

For linux/ubuntu users, you need to install two packages (alternatively use the synaptic package manager)

  sudo apt-get install libopenmpi-dev
  sudo apt-get install openmpi-bin

For OS X users, install brew (after having installed xcode and gcc, needed for the gfortran compiler of openmpi) and then install with brew

   brew install openmpi

When running an executable (code.x), run as

  mpirun -n 10 ./code.x

where we indicate that we want the number of processes to be 10.

Installing MPI and using Qt

With openmpi installed, when using Qt, add to your .pro file the instructions here

You may need to tell Qt where openmpi is stored.

What is Message Passing Interface (MPI)?

MPI is a library, not a language. It specifies the names, calling sequences and results of functions or subroutines to be called from C/C++ or Fortran programs, and the classes and methods that make up the MPI C++ library. The programs that users write in Fortran, C or C++ are compiled with ordinary compilers and linked with the MPI library.

MPI programs should be able to run on all possible machines and run all MPI implementetations without change.

An MPI computation is a collection of processes communicating with messages.

Going Parallel with MPI

Task parallelism: the work of a global problem can be divided into a number of independent tasks, which rarely need to synchronize. Monte Carlo simulations or numerical integration are examples of this.

MPI is a message-passing library where all the routines have corresponding C/C++-binding

   MPI_Command_name

and Fortran-binding (routine names are in uppercase, but can also be in lower case)

   MPI_COMMAND_NAME

MPI is a library

MPI is a library specification for the message passing interface, proposed as a standard.

independent of hardware;
not a language or compiler specification;
not a specific implementation or product.

A message passing standard for portability and ease-of-use. Designed for high performance.

Insert communication and synchronization functions where necessary.

Bindings to MPI routines

MPI is a message-passing library where all the routines have corresponding C/C++-binding

   MPI_Command_name

and Fortran-binding (routine names are in uppercase, but can also be in lower case)

   MPI_COMMAND_NAME

The discussion in these slides focuses on the C++ binding.

Communicator

A group of MPI processes with a name (context).
Any process is identified by its rank. The rank is only meaningful within a particular communicator.
By default the communicator contains all the MPI processes.

  MPI_COMM_WORLD

Mechanism to identify subset of processes.
Promotes modular design of parallel libraries.

Some of the most important MPI functions

$ MPI\_Init $ - initiate an MPI computation
$ MPI\_Finalize $ - terminate the MPI computation and clean up
$ MPI\_Comm\_size $ - how many processes participate in a given MPI communicator?
$ MPI\_Comm\_rank $ - which one am I? (A number between 0 and size-1.)
$ MPI\_Send $ - send a message to a particular process within an MPI communicator
$ MPI\_Recv $ - receive a message from a particular process within an MPI communicator
$ MPI\_reduce $ or $ MPI\_Allreduce $, send and receive messages

The first MPI C/C++ program

Let every process write "Hello world" (oh not this program again!!) on the standard output.

using namespace std;
#include <mpi.h>
#include <iostream>
int main (int nargs, char* args[])
{
int numprocs, my_rank;
//   MPI initializations
MPI_Init (&nargs, &args);
MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
cout << "Hello world, I have  rank " << my_rank << " out of " 
     << numprocs << endl;
//  End MPI
MPI_Finalize ();

The Fortran program

PROGRAM hello
INCLUDE "mpif.h"
INTEGER:: size, my_rank, ierr

CALL  MPI_INIT(ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, my_rank, ierr)
WRITE(*,*)"Hello world, I've rank ",my_rank," out of ",size
CALL MPI_FINALIZE(ierr)

END PROGRAM hello

Note 1

The output to screen is not ordered since all processes are trying to write to screen simultaneously.
It is the operating system which opts for an ordering.
If we wish to have an organized output, starting from the first process, we may rewrite our program as in the next example.

Ordered output with MPIBarrier

int main (int nargs, char* args[])
{
 int numprocs, my_rank, i;
 MPI_Init (&nargs, &args);
 MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
 MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
 for (i = 0; i < numprocs; i++) {}
 MPI_Barrier (MPI_COMM_WORLD);
 if (i == my_rank) {
 cout << "Hello world, I have  rank " << my_rank << 
        " out of " << numprocs << endl;}
      MPI_Finalize ();

Note 2

Here we have used the $ MPI\_Barrier $ function to ensure that that every process has completed its set of instructions in a particular order.
A barrier is a special collective operation that does not allow the processes to continue until all processes in the communicator (here $ MPI\_COMM\_WORLD $) have called $ MPI\_Barrier $.
The barriers make sure that all processes have reached the same point in the code. Many of the collective operations like $ MPI\_ALLREDUCE $ to be discussed later, have the same property; that is, no process can exit the operation until all processes have started.

However, this is slightly more time-consuming since the processes synchronize between themselves as many times as there are processes. In the next Hello world example we use the send and receive functions in order to a have a synchronized action.

Ordered output

.....
int numprocs, my_rank, flag;
MPI_Status status;
MPI_Init (&nargs, &args);
MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
if (my_rank > 0)
MPI_Recv (&flag, 1, MPI_INT, my_rank-1, 100, 
           MPI_COMM_WORLD, &status);
cout << "Hello world, I have  rank " << my_rank << " out of " 
<< numprocs << endl;
if (my_rank < numprocs-1)
MPI_Send (&my_rank, 1, MPI_INT, my_rank+1, 
          100, MPI_COMM_WORLD);
MPI_Finalize ();

Note 3

The basic sending of messages is given by the function $ MPI\_SEND $, which in C/C++ is defined as

int MPI_Send(void *buf, int count, 
             MPI_Datatype datatype, 
             int dest, int tag, MPI_Comm comm)}

This single command allows the passing of any kind of variable, even a large array, to any group of tasks. The variable buf is the variable we wish to send while count is the number of variables we are passing. If we are passing only a single value, this should be 1.

If we transfer an array, it is the overall size of the array. For example, if we want to send a 10 by 10 array, count would be $ 10\times 10=100 $ since we are actually passing 100 values.

Note 4

Once you have sent a message, you must receive it on another task. The function $ MPI\_RECV $ is similar to the send call.

int MPI_Recv( void *buf, int count, MPI_Datatype datatype, 
            int source, 
            int tag, MPI_Comm comm, MPI_Status *status )

The arguments that are different from those in MPI\_SEND are buf which is the name of the variable where you will be storing the received data, source which replaces the destination in the send command. This is the return ID of the sender.

Finally, we have used $ MPI\_Status\_status $, where one can check if the receive was completed.

The output of this code is the same as the previous example, but now process 0 sends a message to process 1, which forwards it further to process 2, and so forth.

Numerical integration in parallel

Integrating $ \pi $

The code example computes $ \pi $ using the trapezoidal rules.
The trapezoidal rule

$$ I=\int_a^bf(x) dx\approx h\left(f(a)/2 + f(a+h) +f(a+2h)+\dots +f(b-h)+ f(b)/2\right). $$

Click on this link for the full program.

Dissection of trapezoidal rule with $ MPI\_reduce $

//    Trapezoidal rule and numerical integration usign MPI
using namespace std;
#include <mpi.h>
#include <iostream>

//     Here we define various functions called by the main program

double int_function(double );
double trapezoidal_rule(double , double , int , double (*)(double));

//   Main function begins here
int main (int nargs, char* args[])
{
  int n, local_n, numprocs, my_rank; 
  double a, b, h, local_a, local_b, total_sum, local_sum;   
  double  time_start, time_end, total_time;

Dissection of trapezoidal rule

  //  MPI initializations
  MPI_Init (&nargs, &args);
  MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
  time_start = MPI_Wtime();
  //  Fixed values for a, b and n 
  a = 0.0 ; b = 1.0;  n = 1000;
  h = (b-a)/n;    // h is the same for all processes 
  local_n = n/numprocs;  
  // make sure n > numprocs, else integer division gives zero
  // Length of each process' interval of
  // integration = local_n*h.  
  local_a = a + my_rank*local_n*h;
  local_b = local_a + local_n*h;

Integrating with MPI

  total_sum = 0.0;
  local_sum = trapezoidal_rule(local_a, local_b, local_n, 
                               &int_function); 
  MPI_Reduce(&local_sum, &total_sum, 1, MPI_DOUBLE, 
              MPI_SUM, 0, MPI_COMM_WORLD);
  time_end = MPI_Wtime();
  total_time = time_end-time_start;
  if ( my_rank == 0) {
    cout << "Trapezoidal rule = " <<  total_sum << endl;
    cout << "Time = " <<  total_time  
         << " on number of processors: "  << numprocs  << endl;
  }
  // End MPI
  MPI_Finalize ();  
  return 0;
}  // end of main program

How do I use $ MPI\_reduce $?

Here we have used

MPI_reduce( void *senddata, void* resultdata, int count, 
     MPI_Datatype datatype, MPI_Op, int root, MPI_Comm comm)

The two variables $ senddata $ and $ resultdata $ are obvious, besides the fact that one sends the address of the variable or the first element of an array. If they are arrays they need to have the same size. The variable $ count $ represents the total dimensionality, 1 in case of just one variable, while $ MPI\_Datatype $ defines the type of variable which is sent and received.

The new feature is $ MPI\_Op $. It defines the type of operation we want to do.

More on $ MPI\_Reduce $

In our case, since we are summing the rectangle contributions from every process we define $ MPI\_Op = MPI\_SUM $. If we have an array or matrix we can search for the largest og smallest element by sending either $ MPI\_MAX $ or $ MPI\_MIN $. If we want the location as well (which array element) we simply transfer $ MPI\_MAXLOC $ or $ MPI\_MINOC $. If we want the product we write $ MPI\_PROD $.

$ MPI\_Allreduce $ is defined as

MPI_Allreduce( void *senddata, void* resultdata, int count, 
          MPI_Datatype datatype, MPI_Op, MPI_Comm comm)

Dissection of trapezoidal rule

We use $ MPI\_reduce $ to collect data from each process. Note also the use of the function $ MPI\_Wtime $.

//  this function defines the function to integrate
double int_function(double x)
{
  double value = 4./(1.+x*x);
  return value;
} // end of function to evaluate

Dissection of trapezoidal rule

//  this function defines the trapezoidal rule
double trapezoidal_rule(double a, double b, int n, 
                         double (*func)(double))
{
  double trapez_sum;
  double fa, fb, x, step;
  int    j;
  step=(b-a)/((double) n);
  fa=(*func)(a)/2. ;
  fb=(*func)(b)/2. ;
  trapez_sum=0.;
  for (j=1; j <= n-1; j++){
    x=j*step+a;
    trapez_sum+=(*func)(x);
  }
  trapez_sum=(trapez_sum+fb+fa)*step;
  return trapez_sum;
}  // end trapezoidal_rule

The quantum dot program for two electrons

// Variational Monte Carlo for atoms with importance sampling, slater det
// Test case for 2-electron quantum dot, no classes using Mersenne-Twister RNG
#include "mpi.h"
#include <cmath>
#include <random>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "vectormatrixclass.h"

using namespace  std;
// output file as global variable
ofstream ofile;  
// the step length and its squared inverse for the second derivative 
//  Here we define global variables  used in various functions
//  These can be changed by using classes
int Dimension = 2; 
int NumberParticles  = 2;  //  we fix also the number of electrons to be 2

// declaration of functions 

// The Mc sampling for the variational Monte Carlo 
void  MonteCarloSampling(int, double &, double &, Vector &);

// The variational wave function
double  WaveFunction(Matrix &, Vector &);

// The local energy 
double  LocalEnergy(Matrix &, Vector &);

// The quantum force
void  QuantumForce(Matrix &, Matrix &, Vector &);


// inline function for single-particle wave function
inline double SPwavefunction(double r, double alpha) { 
   return exp(-alpha*r*0.5);
}

// inline function for derivative of single-particle wave function
inline double DerivativeSPwavefunction(double r, double alpha) { 
  return -r*alpha;
}

// function for absolute value of relative distance
double RelativeDistance(Matrix &r, int i, int j) { 
      double r_ij = 0;  
      for (int k = 0; k < Dimension; k++) { 
	r_ij += (r(i,k)-r(j,k))*(r(i,k)-r(j,k));
      }
      return sqrt(r_ij); 
}

// inline function for derivative of Jastrow factor
inline double JastrowDerivative(Matrix &r, double beta, int i, int j, int k){
  return (r(i,k)-r(j,k))/(RelativeDistance(r, i, j)*pow(1.0+beta*RelativeDistance(r, i, j),2));
}

// function for square of position of single particle
double singleparticle_pos2(Matrix &r, int i) { 
    double r_single_particle = 0;
    for (int j = 0; j < Dimension; j++) { 
      r_single_particle  += r(i,j)*r(i,j);
    }
    return r_single_particle;
}

void lnsrch(int n, Vector &xold, double fold, Vector &g, Vector &p, Vector &x,
		 double *f, double stpmax, int *check, double (*func)(Vector &p));

void dfpmin(Vector &p, int n, double gtol, int *iter, double *fret,
	    double(*func)(Vector &p), void (*dfunc)(Vector &p, Vector &g));

static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)


static double maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
        (maxarg1) : (maxarg2))


// Begin of main program   

int main(int argc, char* argv[])
{

  //  MPI initializations
  int NumberProcesses, MyRank, NumberMCsamples;
  MPI_Init (&argc, &argv);
  MPI_Comm_size (MPI_COMM_WORLD, &NumberProcesses);
  MPI_Comm_rank (MPI_COMM_WORLD, &MyRank);
  double StartTime = MPI_Wtime();
  if (MyRank == 0 && argc <= 1) {
    cout << "Bad Usage: " << argv[0] << 
      " Read also output file on same line and number of Monte Carlo cycles" << endl;
  }
  // Read filename and number of Monte Carlo cycles from the command line
  if (MyRank == 0 && argc > 2) {
    string filename = argv[1]; // first command line argument after name of program
    NumberMCsamples  = atoi(argv[2]);
    string fileout = filename;
    string argument = to_string(NumberMCsamples);
    // Final filename as filename+NumberMCsamples
    fileout.append(argument);
    ofile.open(fileout);
  }
  // broadcast the number of  Monte Carlo samples
  MPI_Bcast (&NumberMCsamples, 1, MPI_INT, 0, MPI_COMM_WORLD);
  // Two variational parameters only
  Vector VariationalParameters(2);
  int TotalNumberMCsamples = NumberMCsamples*NumberProcesses; 
  // Loop over variational parameters
  for (double alpha = 0.5; alpha <= 1.5; alpha +=0.1){
    for (double beta = 0.1; beta <= 0.5; beta +=0.05){
      VariationalParameters(0) = alpha;  // value of alpha
      VariationalParameters(1) = beta;  // value of beta
      //  Do the mc sampling  and accumulate data with MPI_Reduce
      double TotalEnergy, TotalEnergySquared, LocalProcessEnergy, LocalProcessEnergy2;
      LocalProcessEnergy = LocalProcessEnergy2 = 0.0;
      MonteCarloSampling(NumberMCsamples, LocalProcessEnergy, LocalProcessEnergy2, VariationalParameters);
      //  Collect data in total averages
      MPI_Reduce(&LocalProcessEnergy, &TotalEnergy, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
      MPI_Reduce(&LocalProcessEnergy2, &TotalEnergySquared, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
      // Print out results  in case of Master node, set to MyRank = 0
      if ( MyRank == 0) {
	double Energy = TotalEnergy/( (double)NumberProcesses);
	double Variance = TotalEnergySquared/( (double)NumberProcesses)-Energy*Energy;
	double StandardDeviation = sqrt(Variance/((double)TotalNumberMCsamples)); // over optimistic error
	ofile << setiosflags(ios::showpoint | ios::uppercase);
	ofile << setw(15) << setprecision(8) << VariationalParameters(0);
	ofile << setw(15) << setprecision(8) << VariationalParameters(1);
	ofile << setw(15) << setprecision(8) << Energy;
	ofile << setw(15) << setprecision(8) << Variance;
	ofile << setw(15) << setprecision(8) << StandardDeviation << endl;
      }
    }
  }
  double EndTime = MPI_Wtime();
  double TotalTime = EndTime-StartTime;
  if ( MyRank == 0 )  cout << "Time = " <<  TotalTime  << " on number of processors: "  << NumberProcesses  << endl;
  if (MyRank == 0)  ofile.close();  // close output file
  // End MPI
  MPI_Finalize ();  
  return 0;
}  //  end of main function


// Monte Carlo sampling with the Metropolis algorithm  

void MonteCarloSampling(int NumberMCsamples, double &cumulative_e, double &cumulative_e2, Vector &VariationalParameters)
{

 // Initialize the seed and call the Mersienne algo
  std::random_device rd;
  std::mt19937_64 gen(rd());
  // Set up the uniform distribution for x \in [[0, 1]
  std::uniform_real_distribution<double> UniformNumberGenerator(0.0,1.0);
  std::normal_distribution<double> Normaldistribution(0.0,1.0);
  // diffusion constant from Schroedinger equation
  double D = 0.5; 
  double timestep = 0.05;  //  we fix the time step  for the gaussian deviate
  // allocate matrices which contain the position of the particles  
  Matrix OldPosition( NumberParticles, Dimension), NewPosition( NumberParticles, Dimension);
  Matrix OldQuantumForce(NumberParticles, Dimension), NewQuantumForce(NumberParticles, Dimension);
  double Energy = 0.0; double EnergySquared = 0.0; double DeltaE = 0.0;
  //  initial trial positions
  for (int i = 0; i < NumberParticles; i++) { 
    for (int j = 0; j < Dimension; j++) {
      OldPosition(i,j) = Normaldistribution(gen)*sqrt(timestep);
    }
  }
  double OldWaveFunction = WaveFunction(OldPosition, VariationalParameters);
  QuantumForce(OldPosition, OldQuantumForce, VariationalParameters);
  // loop over monte carlo cycles 
  for (int cycles = 1; cycles <= NumberMCsamples; cycles++){ 
    // new position 
    for (int i = 0; i < NumberParticles; i++) { 
      for (int j = 0; j < Dimension; j++) {
	// gaussian deviate to compute new positions using a given timestep
	NewPosition(i,j) = OldPosition(i,j) + Normaldistribution(gen)*sqrt(timestep)+OldQuantumForce(i,j)*timestep*D;
	//	NewPosition(i,j) = OldPosition(i,j) + gaussian_deviate(&idum)*sqrt(timestep)+OldQuantumForce(i,j)*timestep*D;
      }  
      //  for the other particles we need to set the position to the old position since
      //  we move only one particle at the time
      for (int k = 0; k < NumberParticles; k++) {
	if ( k != i) {
	  for (int j = 0; j < Dimension; j++) {
	    NewPosition(k,j) = OldPosition(k,j);
	  }
	} 
      }
      double NewWaveFunction = WaveFunction(NewPosition, VariationalParameters); 
      QuantumForce(NewPosition, NewQuantumForce, VariationalParameters);
      //  we compute the log of the ratio of the greens functions to be used in the 
      //  Metropolis-Hastings algorithm
      double GreensFunction = 0.0;            
      for (int j = 0; j < Dimension; j++) {
	GreensFunction += 0.5*(OldQuantumForce(i,j)+NewQuantumForce(i,j))*
	  (D*timestep*0.5*(OldQuantumForce(i,j)-NewQuantumForce(i,j))-NewPosition(i,j)+OldPosition(i,j));
      }
      GreensFunction = exp(GreensFunction);
      // The Metropolis test is performed by moving one particle at the time
      if(UniformNumberGenerator(gen) <= GreensFunction*NewWaveFunction*NewWaveFunction/OldWaveFunction/OldWaveFunction ) { 
	for (int  j = 0; j < Dimension; j++) {
	  OldPosition(i,j) = NewPosition(i,j);
	  OldQuantumForce(i,j) = NewQuantumForce(i,j);
	}
	OldWaveFunction = NewWaveFunction;
      }
    }  //  end of loop over particles
    // compute local energy  
    double DeltaE = LocalEnergy(OldPosition, VariationalParameters);
    // update energies
    Energy += DeltaE;
    EnergySquared += DeltaE*DeltaE;
  }   // end of loop over MC trials   
  // update the energy average and its squared 
  cumulative_e = Energy/NumberMCsamples;
  cumulative_e2 = EnergySquared/NumberMCsamples;
}   // end MonteCarloSampling function  


// Function to compute the squared wave function and the quantum force

double  WaveFunction(Matrix &r, Vector &VariationalParameters)
{
  double wf = 0.0;
  // full Slater determinant for two particles, replace with Slater det for more particles 
  wf  = SPwavefunction(singleparticle_pos2(r, 0), VariationalParameters(0))*SPwavefunction(singleparticle_pos2(r, 1),VariationalParameters(0));
  // contribution from Jastrow factor
  for (int i = 0; i < NumberParticles-1; i++) { 
    for (int j = i+1; j < NumberParticles; j++) {
      wf *= exp(RelativeDistance(r, i, j)/((1.0+VariationalParameters(1)*RelativeDistance(r, i, j))));
    }
  }
  return wf;
}

// Function to calculate the local energy without numerical derivation of kinetic energy

double  LocalEnergy(Matrix &r, Vector &VariationalParameters)
{

  // compute the kinetic and potential energy from the single-particle part
  // for a many-electron system this has to be replaced by a Slater determinant
  // The absolute value of the interparticle length
  Matrix length( NumberParticles, NumberParticles);
  // Set up interparticle distance
  for (int i = 0; i < NumberParticles-1; i++) { 
    for(int j = i+1; j < NumberParticles; j++){
      length(i,j) = RelativeDistance(r, i, j);
      length(j,i) =  length(i,j);
    }
  }
  double KineticEnergy = 0.0;
  // Set up kinetic energy from Slater and Jastrow terms
  for (int i = 0; i < NumberParticles; i++) { 
    for (int k = 0; k < Dimension; k++) {
      double sum1 = 0.0; 
      for(int j = 0; j < NumberParticles; j++){
	if ( j != i) {
	  sum1 += JastrowDerivative(r, VariationalParameters(1), i, j, k);
	}
      }
      KineticEnergy += (sum1+DerivativeSPwavefunction(r(i,k),VariationalParameters(0)))*(sum1+DerivativeSPwavefunction(r(i,k),VariationalParameters(0)));
    }
  }
  KineticEnergy += -2*VariationalParameters(0)*NumberParticles;
  for (int i = 0; i < NumberParticles-1; i++) {
      for (int j = i+1; j < NumberParticles; j++) {
        KineticEnergy += 2.0/(pow(1.0 + VariationalParameters(1)*length(i,j),2))*(1.0/length(i,j)-2*VariationalParameters(1)/(1+VariationalParameters(1)*length(i,j)) );
      }
  }
  KineticEnergy *= -0.5;
  // Set up potential energy, external potential + eventual electron-electron repulsion
  double PotentialEnergy = 0;
  for (int i = 0; i < NumberParticles; i++) { 
    double DistanceSquared = singleparticle_pos2(r, i);
    PotentialEnergy += 0.5*DistanceSquared;  // sp energy HO part, note it has the oscillator frequency set to 1!
  }
  // Add the electron-electron repulsion
  for (int i = 0; i < NumberParticles-1; i++) { 
    for (int j = i+1; j < NumberParticles; j++) {
      PotentialEnergy += 1.0/length(i,j);          
    }
  }
  double LocalE = KineticEnergy+PotentialEnergy;
  return LocalE;
}

// Compute the analytical expression for the quantum force
void  QuantumForce(Matrix &r, Matrix &qforce, Vector &VariationalParameters)
{
  // compute the first derivative 
  for (int i = 0; i < NumberParticles; i++) {
    for (int k = 0; k < Dimension; k++) {
      // single-particle part, replace with Slater det for larger systems
      double sppart = DerivativeSPwavefunction(r(i,k),VariationalParameters(0));
      //  Jastrow factor contribution
      double Jsum = 0.0;
      for (int j = 0; j < NumberParticles; j++) {
	if ( j != i) {
	  Jsum += JastrowDerivative(r, VariationalParameters(1), i, j, k);
	}
      }
      qforce(i,k) = 2.0*(Jsum+sppart);
    }
  }
} // end of QuantumForce function


#define ITMAX 200
#define EPS 3.0e-8
#define TOLX (4*EPS)
#define STPMX 100.0

void dfpmin(Vector &p, int n, double gtol, int *iter, double *fret,
	    double(*func)(Vector &p), void (*dfunc)(Vector &p, Vector &g))
{

  int check,i,its,j;
  double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
  Vector dg(n), g(n), hdg(n), pnew(n), xi(n);
  Matrix hessian(n,n);

  fp=(*func)(p);
  (*dfunc)(p,g);
  for (i = 0;i < n;i++) {
    for (j = 0; j< n;j++) hessian(i,j)=0.0;
    hessian(i,i)=1.0;
    xi(i) = -g(i);
    sum += p(i)*p(i);
  }
  stpmax=STPMX*FMAX(sqrt(sum),(double)n);
  for (its=1;its<=ITMAX;its++) {
    *iter=its;
    lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,&check,func);
    fp = *fret;
    for (i = 0; i< n;i++) {
      xi(i)=pnew(i)-p(i);
      p(i)=pnew(i);
    }
    test=0.0;
    for (i = 0;i< n;i++) {
      temp=fabs(xi(i))/FMAX(fabs(p(i)),1.0);
      if (temp > test) test=temp;
    }
    if (test < TOLX) {
      return;
    }
    for (i=0;i<n;i++) dg(i)=g(i);
    (*dfunc)(p,g);
    test=0.0;
    den=FMAX(*fret,1.0);
    for (i=0;i<n;i++) {
      temp=fabs(g(i))*FMAX(fabs(p(i)),1.0)/den;
      if (temp > test) test=temp;
    }
    if (test < gtol) {
      return;
    }
    for (i=0;i<n;i++) dg(i)=g(i)-dg(i);
    for (i=0;i<n;i++) {
      hdg(i)=0.0;
      for (j=0;j<n;j++) hdg(i) += hessian(i,j)*dg(j);
    }
    fac=fae=sumdg=sumxi=0.0;
    for (i=0;i<n;i++) {
      fac += dg(i)*xi(i);
      fae += dg(i)*hdg(i);
      sumdg += SQR(dg(i));
      sumxi += SQR(xi(i));
    }
    if (fac*fac > EPS*sumdg*sumxi) {
      fac=1.0/fac;
      fad=1.0/fae;
      for (i=0;i<n;i++) dg(i)=fac*xi(i)-fad*hdg(i);
      for (i=0;i<n;i++) {
	for (j=0;j<n;j++) {
	  hessian(i,j) += fac*xi(i)*xi(j)
	    -fad*hdg(i)*hdg(j)+fae*dg(i)*dg(j);
	}
      }
    }
    for (i=0;i<n;i++) {
      xi(i)=0.0;
      for (j=0;j<n;j++) xi(i) -= hessian(i,j)*g(j);
    }
  }
  cout << "too many iterations in dfpmin" << endl;
}
#undef ITMAX
#undef EPS
#undef TOLX
#undef STPMX

#define ALF 1.0e-4
#define TOLX 1.0e-7

void lnsrch(int n, Vector &xold, double fold, Vector &g, Vector &p, Vector &x,
	    double *f, double stpmax, int *check, double (*func)(Vector &p))
{
  int i;
  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,
    test,tmplam;

  *check=0;
  for (sum=0.0,i=0;i<n;i++) sum += p(i)*p(i);
  sum=sqrt(sum);
  if (sum > stpmax)
    for (i=0;i<n;i++) p(i) *= stpmax/sum;
  for (slope=0.0,i=0;i<n;i++)
    slope += g(i)*p(i);
  test=0.0;
  for (i=0;i<n;i++) {
    temp=fabs(p(i))/FMAX(fabs(xold(i)),1.0);
    if (temp > test) test=temp;
  }
  alamin=TOLX/test;
  alam=1.0;
  for (;;) {
    for (i=0;i<n;i++) x(i)=xold(i)+alam*p(i);
    *f=(*func)(x);
    if (alam < alamin) {
      for (i=0;i<n;i++) x(i)=xold(i);
      *check=1;
      return;
    } else if (*f <= fold+ALF*alam*slope) return;
    else {
      if (alam == 1.0)
	tmplam = -slope/(2.0*(*f-fold-slope));
      else {
	rhs1 = *f-fold-alam*slope;
	rhs2=f2-fold2-alam2*slope;
	a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
	b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
	if (a == 0.0) tmplam = -slope/(2.0*b);
	else {
	  disc=b*b-3.0*a*slope;
	  if (disc<0.0) cout << "Roundoff problem in lnsrch." << endl;
	  else tmplam=(-b+sqrt(disc))/(3.0*a);
	}
	if (tmplam>0.5*alam)
	  tmplam=0.5*alam;
      }
    }
    alam2=alam;
    f2 = *f;
    fold2=fold;
    alam=FMAX(tmplam,0.1*alam);
  }
}
#undef ALF
#undef TOLX

What is OpenMP

OpenMP provides high-level thread programming
Multiple cooperating threads are allowed to run simultaneously
Threads are created and destroyed dynamically in a fork-join pattern

An OpenMP program consists of a number of parallel regions
Between two parallel regions there is only one master thread
In the beginning of a parallel region, a team of new threads is spawned

The newly spawned threads work simultaneously with the master thread
At the end of a parallel region, the new threads are destroyed

Many good tutorials online and excellent textbook

Using OpenMP, by B. Chapman, G. Jost, and A. van der Pas
Many tutorials online like OpenMP official site

Getting started, things to remember

Remember the header file

#include <omp.h>

Insert compiler directives in C++ syntax as

#pragma omp...

Compile with for example c++ -fopenmp code.cpp
Execute

Remember to assign the environment variable OMP NUM THREADS
It specifies the total number of threads inside a parallel region, if not otherwise overwritten

OpenMP syntax

Different OpenMP styles of parallelism

General code structure

#include <omp.h>
main ()
{
int var1, var2, var3;
/* serial code */
/* ... */
/* start of a parallel region */
#pragma omp parallel private(var1, var2) shared(var3)
{
/* ... */
}
/* more serial code */
/* ... */
/* another parallel region */
#pragma omp parallel
{
/* ... */
}
}

Parallel region

A parallel region is a block of code that is executed by a team of threads
The following compiler directive creates a parallel region

#pragma omp parallel { ... }

Clauses can be added at the end of the directive
Most often used clauses:

default(shared) or default(none)
public(list of variables)
private(list of variables)

Hello world, not again, please!

#include <omp.h>
#include <cstdio>
int main (int argc, char *argv[])
{
int th_id, nthreads;
#pragma omp parallel private(th_id) shared(nthreads)
{
th_id = omp_get_thread_num();
printf("Hello World from thread %d\n", th_id);
#pragma omp barrier
if ( th_id == 0 ) {
nthreads = omp_get_num_threads();
printf("There are %d threads\n",nthreads);
}
}
return 0;
}

Hello world, yet another variant

#include <cstdio>
#include <omp.h>
int main(int argc, char *argv[]) 
{
 omp_set_num_threads(4); 
#pragma omp parallel
 {
   int id = omp_get_thread_num();
   int nproc = omp_get_num_threads(); 
   cout << "Hello world with id number and processes " <<  id <<  nproc << endl;
 } 
return 0;
}

Variables declared outside of the parallel region are shared by all threads If a variable like id is declared outside of the

#pragma omp parallel,

it would have been shared by various the threads, possibly causing erroneous output

Why? What would go wrong? Why do we add possibly?

Important OpenMP library routines

int omp get num threads (), returns the number of threads inside a parallel region
int omp get thread num (), returns the a thread for each thread inside a parallel region
void omp set num threads (int), sets the number of threads to be used
void omp set nested (int), turns nested parallelism on/off

Private variables

Private clause can be used to make thread- private versions of such variables:

#pragma omp parallel private(id)
{
 int id = omp_get_thread_num();
 cout << "My thread num" << id << endl; 
}

What is their value on entry? Exit?
OpenMP provides ways to control that
Can use default(none) to require the sharing of each variable to be described

Master region

It is often useful to have only one thread execute some of the code in a parallel region. I/O statements are a common example

#pragma omp parallel 
{
  #pragma omp master
   {
      int id = omp_get_thread_num();
      cout << "My thread num" << id << endl; 
   } 
}

Parallel for loop

Inside a parallel region, the following compiler directive can be used to parallelize a for-loop:

#pragma omp for

Clauses can be added, such as

schedule(static, chunk size)
schedule(dynamic, chunk size)
schedule(guided, chunk size) (non-deterministic allocation)
schedule(runtime)
private(list of variables)
reduction(operator:variable)
nowait

Parallel computations and loops

OpenMP provides an easy way to parallelize a loop

#pragma omp parallel for
  for (i=0; i<n; i++) c[i] = a[i];

OpenMP handles index variable (no need to declare in for loop or make private)

Which thread does which values? Several options.

Scheduling of loop computations

We can let the OpenMP runtime decide. The decision is about how the loop iterates are scheduled and OpenMP defines three choices of loop scheduling:

Static: Predefined at compile time. Lowest overhead, predictable
Dynamic: Selection made at runtime
Guided: Special case of dynamic; attempts to reduce overhead

Example code for loop scheduling

#include <omp.h>
#define CHUNKSIZE 100
#define N 1000
int main (int argc, char *argv[])
{
int i, chunk;
float a[N], b[N], c[N];
for (i=0; i < N; i++) a[i] = b[i] = i * 1.0;
chunk = CHUNKSIZE;
#pragma omp parallel shared(a,b,c,chunk) private(i)
{
#pragma omp for schedule(dynamic,chunk)
for (i=0; i < N; i++) c[i] = a[i] + b[i];
} /* end of parallel region */
}

Example code for loop scheduling, guided instead of dynamic

#include <omp.h>
#define CHUNKSIZE 100
#define N 1000
int main (int argc, char *argv[])
{
int i, chunk;
float a[N], b[N], c[N];
for (i=0; i < N; i++) a[i] = b[i] = i * 1.0;
chunk = CHUNKSIZE;
#pragma omp parallel shared(a,b,c,chunk) private(i)
{
#pragma omp for schedule(guided,chunk)
for (i=0; i < N; i++) c[i] = a[i] + b[i];
} /* end of parallel region */
}

Week 11, March 11-15: Resampling Techniques, Bootstrap and Blocking

March 11-15

Overview of week 11, March 11-15

Why resampling methods ?

Statistical analysis

And why do we use such methods?

Central limit theorem

Further remarks

Running many measurements

Adding more definitions

Further rewriting

The covariance term

Rewriting the covariance term

Introducing the correlation function

Computing the correlation function

Resampling methods: Blocking

Why blocking?

Blocking Transformations

Blocking transformations

Blocking Transformations

Blocking Transformations, getting there

Blocking Transformations, final expressions

More on the blocking method

Example code form last week

Resampling analysis

Content

Optimization and profiling

More on optimization

Optimization and profiling

Optimization and debugging

Other hints

Vectorization and the basic idea behind parallel computing

A rough classification of hardware models

Shared memory and distributed memory

Different parallel programming paradigms

Different parallel programming paradigms

What is vectorization?

Number of elements that can acted upon

Number of elements that can acted upon, examples

Operation counts for scalar operation

Number of elements that can acted upon, examples

Number of operations when vectorized

Compiling with and without vectorization

Compiling with and without vectorization using clang

Automatic vectorization and vectorization inhibitors, criteria

Automatic vectorization and vectorization inhibitors, exit criteria

Automatic vectorization and vectorization inhibitors, straight-line code

Automatic vectorization and vectorization inhibitors, nested loops

Automatic vectorization and vectorization inhibitors, function calls

Automatic vectorization and vectorization inhibitors, data dependencies

Automatic vectorization and vectorization inhibitors, more data dependencies

Automatic vectorization and vectorization inhibitors, memory stride

Memory management

Memory and communication

Measuring performance

Problems with measuring time

Problems with cold start

Problems with smart compilers

Problems with interference

Problems with measuring performance

Thomas algorithm for tridiagonal linear algebra equations

Thomas algorithm, forward substitution

Thomas algorithm, backward substitution

Thomas algorithm and counting of operations (floating point and memory)

How do we define speedup? Simplest form

How do we define speedup? Correct baseline

Parallel speedup

Speedup and memory

Upper bounds on speedup

Amdahl's law

How much is parallelizable

Today's situation of parallel computing

Overhead present in parallel computing

Parallelizing a sequential algorithm

Strategies

How do I run MPI on a PC/Laptop? MPI

Can I do it on my own PC/laptop? OpenMP installation

Installing MPI

Installing MPI and using Qt

What is Message Passing Interface (MPI)?