Introducción a la Computación Científica y de Alto Rendimiento

• https://tinyurl.com/IntroSciCompHPC

• 2023-I: https://www.youtube.com/playlist?list=PLxbXsIEaf05rKMvhEhoRLHL-yJypwsJPK
• 2022-I: https://www.youtube.com/playlist?list=PLxbXsIEaf05oC3xW7BpAPC-TECWrqV2Ni

• C++ programming : https://www.youtube.com/playlist?list=PLxbXsIEaf05qKMpNR69p4iUerDL9mvcE4

• https://swcarpentry.github.io/shell-novice/
• https://github.com/phyver/GameShell
• https://overthewire.org/wargames/bandit/

• https://swcarpentry.github.io/git-novice/
• http://swcarpentry.github.io/git-novice-es
• https://ohmygit.org/

• We use g++, the gnu c++ compiler
• The most basic form to use it is :

     g++ filename.cpp

If the compilation is successful, this produces a file called a.out on the current directory. You can execute it as

   ./a.out

where ./ means “the current directory”. If you want to change the name of the executable, you can compile as

   g++ filename.cpp -o executablename.x

Replace executablename.x with the name you want for your executable.

 1: 
 2:      // imprima los numeros del 1 al 10 suando while
 3: 
 4:      #include <iostream> 
 5:      using namespace std;
 6: 
 7:      int main(void) 
 8:      {
 9:        int n;
10:        
11:        n = 1;
12:        while (n <= 10) {
13:          cout << n << endl;
14:          ++n;  
15:        } 
16:        
17:        return 0;
18:      }

#+RESULTS

We will add 1 until some error occurs

#include <cstdio>

int main(void)
{
  int a = 1;
  while(a > 0) {
    a *= 2 ;
    std::printf("%10d\n", a);
  }
}

Is there any computation difference between the following sums?

Implement and plot the relative difference among them as a function of \(N\).

(Check presentation for a practical example)

• Use numerical libraries: lapack, eigen, standard c++ lib (https://en.cppreference.com/w/, https://hackingcpp.com/cpp/cheat_sheets.html)
• Analize extreme cases, try to rewrite the expressions: http://herbie.uwplse.org/
• Minimize the use of substractions of similar numbers, operations between very large and/or very small numbers, normalize your models
• Some interesting refs:
  • https://phys.org/news/2021-07-rounding-errors-stopwatches-wrong-winners.html
  • https://docs.oracle.com/cd/E19957-01/800-7895/800-7895.pdf

Note: For more information, please see:

• Gnu make tutorial
• Software Carpentry Tutorial .
• https://kbroman.org/minimal_make/
• https://swcarpentry.github.io/make-novice/
• First linux version: https://elixir.bootlin.com/linux/0.01/source/kernel/Makefile (more info at https://lwn.net/SubscriberLink/928581/841b747332791ac4/)

A makefile allows you to atomatize the compilation process by following given rules, and for large projects this could actually speed up significantly the compilation process. You can also use a Makefile for latex projects and other related task where you want to reload and create automatically an updated version of a document when you have updated a small part. A rule in a makefile is written as

target : dependencies separated by space
    rule to build target
#sumupdown.x : sumupdown.cpp
#   g++ sumupdown.cpp -o sumupdown.x

It is important to note that the beginning of the rule is indented with real tab, not with spaces.

Imagine that you want to automatize the creation of the figure of the sumup and sumdown exercise from the numerical errors section. In this case you should think about the dependencies: the latest figure, sumupdown.pdf will depend on the a given gnuplot (or python) script (plot.gp) that read some data from a given file (data.txt), so any change on those files should trigger and update on the figure. Furhtermore, the data file will actually depend on the source (sumupdown.cpp) and its related executable (sumupdown.x), since any change there may imply new data. And actually the executable depends directly on the source code. We could express this dependencies as

fig.pdf: script.gp data.txt
    gnuplot script.gp

data.txt: sumupdown.x
    ./sumupdown.x > data.txt

sumupdown.x: sumupdown.cpp
    g++ sumupdown.cpp -o sumupdown.x

sumupdown.pdf: datos.txt script.gp
    gnuplot script.gp

The actual order is not important. Every time you run make it will check if any target is outdated and run the command needed. The follwing code shows a simple solution for the sum up and down

#include <iostream>
#include <cmath>

float sumup(int nterms);
float sumdown(int nterms);

int main(int argc, char **argv)
{
  std::cout.precision(6); std::cout.setf(std::ios::scientific);
  for (int ii = 1; ii <= 10000; ++ii) {
    double sum1 = sumup(ii);
    double sum2 = sumdown(ii);
    std::cout << ii << "\t" << sum1 << "\t" << sum2
              << "\t" << std::fabs(sum1-sum2)/sum2 << "\n";
  }
  
  return 0;
}

float sumup(int nterms)
{
  float result = 0;
  for(int n = 1; n <= nterms; n++) {
    result += 1.0/n;
  }
  return result;
}
float sumdown(int nterms)
{
  float result = 0;
  for(int n = nterms; n >= 1; n--) {
    result += 1.0/n;
  }
  return result;
}

And this is a simple gnuplot script example

set term pdf
set out "sumupdown.pdf"
# set xlabel "Nterms"
# set ylabel "Percentual difference"
plot 'data.txt' u 1:4 w lp pt 4

Now just run make and see what happens. Then uncomment, for example, the commented lines inside the gnuplot script and run make again. Check that make just reruns the needed rule, it does not recompile the whole program.

Actually, the makefile can be simplified and generalized using variables

fig.pdf: script.gp data.txt
    gnuplot $^

data.txt: sumupdown.x
    ./$^ > $@

%.x: %.cpp
    g++ $^ -o $@

• just : https://github.com/casey/just, https://cheatography.com/linux-china/cheat-sheets/justfile/
• cmake: https://cmake.org/, https://crascit.com/professional-cmake/
• xmake: https://xmake.io/#/
• ninja, bazel, ...

For this example let’s assume that you have two functions, foo and bar, that you want to use repeatedly on several projects so you have moved them to their corresponding header and implementations files as

// file foo.h
#include <iostream>
void foo(void);

  // file foo.cpp
  #include "foo.h"

  void foo(void){
    std::cout << "Inside foo\n";
  }

// file bar.h
#include <iostream>
void bar(void);

  // file bar.cpp
  #include "bar.h"

  void bar(void){
    std::cout << "Inside bar\n";
  }

and the main file,

  #include "foo.h"
  #include "bar.h"

  int main(void)
  {
    foo();
    bar();

    return 0;
  }

If you want to compile these files, you need to run something like

# this creates the object foo.o
g++ -c foo.cpp
# this creates the object bar.o
g++ -c bar.cpp
# this links with main and creates executable main.x
g++ main.cpp foo.o bar.o -o main.x

Of course, all this process can be automatized in a single shell script but, if you change only foo.cpp, the shell script will run again all the commands, so many non needed compilations are performed. A makefile optimizes this process. Let’s write our first makefile for the problem at hand:

main.x : foo.o bar.o
  g++ -o main.x foo.o bar.o main.cpp

Here you see a special structure : a target name (main.x), a pre-requisites list (~foo.o and bar.o), and the rule to make the target. The rule MUST be indented with real tab. Save the file as makefile-v1 and run the command

make -f makefile-v1

and you will see that make will create the objects and the executable. This is great! Even more, the objects are created automatically. These is done through something called automatic rules.

Now change a little bit the file foo.cpp, add some new printing, and re-run make as before. What do you see? Now make is recompiling only the object which has changed, not all of them, so the re-compilation is much faster!.

After seeing this advantages, let’s try to generalize and use more useful syntax for make. Fir instance, you can specify variables, which can be overridden on the command line. Something like

CXX=g++
CXXFLAGS=-I.

main.x: foo.o bar.o
  $(CXX) $(CXXFLAGS) -o main.x main.cpp foo.o bar.o

Compile again and test. Everything should work.

Now, we will specify a generic rule to create the .o objects, and also the name of the objects in a variable. Furthermore, we use $@ to specify the name of the target, and $^ which symbolizes all the items in the dependency list

CXX = g++
CXXFLAGS = -I.
OBJ = foo.o bar.o

main.x: $(OBJ)
  $(CXX) $(CXXFLAGS) -o $@ $^ main.cpp

Save and run it, you will still get the same but the makefile is becoming more generic, and therefore more useful.

Now let’s specify the generic rule to create the .o files, and also specify that if we change the .h files then a recompilation is needed (for now our makefile only detect changes in .cpp files):

CXX = g++
CXXFLAGS = -I.
OBJ = foo.o bar.o
DEPS = foo.h bar.h

%.o: %.cpp $(DEPS)
  $(CXX) -c -o $@ $< $(CXXFLAGS)

main.x: $(OBJ)
  $(CXX) $(CXXFLAGS) -o $@ $^ main.cpp

#+RESULTS Where we have specify the dependency on DEPS, and also we are using a generic rule to tell how to create any .o file (exemplified with the generic %.o symbol), from any corresponding %.c file, and $< means the first item on the deps list. Again, save it and run it.

You can rules to clean the directory (to erase object files, for example), to add libs, to put the headers on another directory, etc. For example, the following adds a phony rule (which does not create any target) to clean the directory

CXX = g++
CXXFLAGS = -I.
OBJ = foo.o bar.o
DEPS = foo.h bar.h

%.o: %.cpp $(DEPS)
  $(CXX) -c -o $@ $< $(CXXFLAGS)

main.x: $(OBJ)
  $(CXX) $(CXXFLAGS) -o $@ $^ main.cpp

.PHONY: clean
clean:
  rm -f *.o *~ *.x

In order to run this rule, you need to write

make -f makefile-v5 clean

And, finally, a more complete example with comments and new syntax, which you can adapt for your own needs

CXX = g++
CXXFLAGS = -I.
LDFLAGS =
SOURCES = main.cpp foo.cpp bar.cpp
OBJ = $(SOURCES:.cpp=.o) # extracts automatically the objects names
DEPS = foo.h bar.h

all : main.x $(SOURCES) $(DEPS) # this is the default target

main.x: $(OBJ)
  @echo "Creating main executable ..."
  $(CXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS)

.cpp.o:
  $(CXX) -c -o $@ $< $(CXXFLAGS)

.PHONY: clean
clean:
  rm -f *.o *~ *.x

1. Take the codes for overflow and underflow and automatize the detection of overflow with the standard function isinf.
2. When they are detecting overflow automatically, stopping and printing where the overflow occurs, put them in separated files (.h and .cpp).
3. Then create a main file which calls those functions.
4. Now write a makefile to automatize the compilation. Use the latest version of the example Makefile.
5. When everything is working send it to assigned repository.
6. Do the same now including the code which computes the machine eps, using its own .h and .cpp .
7. If you use latex, how will you write a Makefile for that? Well, at the end you can use latexmk, of course. Try to think on other uses for make (like automatic zipping and mailing of some important file, )

It is important for you to know a little how your operating system and the shell look for commands. In general, the PATH variable stores the directories where the shell interpreter will look for a given command. To check its contents, you can run the following,

   echo $PATH

If, for instance, you want to add another directory to the PATH, like the directory $HOME/local/bin, you can run the following

   export PATH=$PATH:$HOME/bin

This appends the special directory to the old content of the PATH.

When you want to compile a given program that uses some other libraries, you must specify any extra folder to look for include files, done with the -I flag, and for libraries, done with the -L and -l flags. For instance, let’s assume that you installed some programs in your HOME/local, the include files inside $HOME/local/include, and the libraries inside $HOME/local/lib. If you want to tell the compiler to use those, you must compile as

   g++ -I $HOME/local/include -L $HOME/local/lib filename.cpp -llibname

Finally, whenever you are installing a program from source you must be aware that this reduces to basically the following steps:

1. Download and untar the file. Enter the unpacked directory.
2. Read the README/INSTALL files for any important info.

3. If the program uses cmake, create a build dir and then use cmake to generate the Makefiles:

         mkdir build
         cd build
         cmake ../ -DCMAKE_INSTALL_PREFIX=$HOME/local
   

   On the other hand, if the program uses configure, then configure the system to install on the required path

         ./configure --prefix=$HOME/local
   

4. Compile and install the program, maybe using many threads

         make -j 4 # uses for threads to compile
         make install
   

   Done. The program is installed.

   Finally, when compiling, do not forget to use the flags -L and -I appropriately.

Setting all the flags and making sure to use the right version is sometimes difficult, so tools like spack aim to manage and simplify this.

If you are used to apt-get or something related, you can use the package manager to check. But, in general, you can check the versions by looking at the appropriate places on your system. Typically, if you are looking for a library, they are installed under /usr/lib or /usr/local/lib, while include files are installed under /usr/include or /usr/local/include . For instance, if you are looking for library foo, then you should look for the file libfoo.a or libfoo.so . One useful utility for this is the command locate or find .

  locate libfoo
  find /usr/lib -iname "*libfoo*"

Execute these commands to check the actual versions for fftw, eigen, and git. What versions do you have? If you a re looking for a program, or an specific version of program, you must check if the program exists by executing it. For command line programs you usually can check the version by using the following

   programname --version

where programname is the name of the command.

In this case we will install everything under the $HOME/local subdirectory inside your home, so please create it. Remember that the symbol $HOME means your home directory. The utilies will then create the appropriate folders there. NOTE: Better use the $HOME var instead of ~, which is another way to state the name of your home.

1. Download the source from the project page. This normally implies downloading a tarball (file ending with .tar.gz or .tar.bz2) .

2. Un-compress the downloaded file. For a tarball, the command will be

         tar xf filename.tar.gz
   

3. Enter to the newly uncompressed folder (almost always usually cd filename).
4. READ the README and/or the INSTALL file to check for important info regarding the program. SOmetimes these files tell you that the installation is special and that you are required to follow some special steps (that will happen with voro++ )

5. CONFIGURATION: You have two options, each one independent on the other:

   1. If the program has a configure script, then just run

        ./configure --help
   

   to check all the available options. Since we want to install on the $HOME/local directory, then we need to run

        ./configure --prefix=$HOME/local
   

   If you don’t specify the prefix, then the program will be installed on the /usr/bin or /usr/local/bin directories, whatever is the default. If these commands ends successfully, it will print some info to the screen and will tell you to go to the next step. Otherwise you will need to read the log and fix the errors (like installing a dependency).

   1. If the program uses cmake, a makefile generator and

   configurator, then you need to do the following:

        mkdir build # creates a build directory to put there the temporary built files
        cd build 
        cmake ../ -DCMAKE_INSTALL_PREFIX:PATH=$HOME/local # configure the building process for the source code located on the parent directory
   

6. COMPILATION: Now that you have configured your installation, you need to compile by using the GNU make utility (Note: All this build utilities come from the gnu organization and are free software as in freedom). If you have several cores, you can use them in parallel, assuming the that the Makefile and your make versions supports it:

         make -j 3 # for three cores, but, if you are unsure, just use one core.
   

   Any errors in this stage should be fixed before going to the next one.

7. INSTALLATION After successful compilation, you can install by using

         make install
   

   This will install the program (libraries, binaries, include files, manual files, etc) onto the prefix directory. If you want to instll system-wide (you did not set the prefix), then you need to use sudo make install . In this case you don’t need sudo since you are installing on your own home.

8. TESTING In this case use a program to test your installation. When you compile your program and you want to use the version that you installed, you need to tell the compiler where to find the libraries/includes, so you need something like

         g++ -L $HOME/local/lib -I $HOME/local/include  programname.cpp -llibname
   

   • -L $HOME/local/lib tells the compiler to look for libraries on

   the $HOME/local/lib directory.

   • -I $HOME/local/include tells the compiler to look for include

   files on the $HOME/local/include directory.

   • -llibname tells the compiler to link with the given

   library. Sometimes is not needed. Sometimes is crucial. Be careful, if your library is called libfftw, you need to write -lfftw, not -llibfftw.

For each of the proposed utilities written at the beginning, follow the next procedure:

• Check the installed version number and compare with the latest one.
• Install the latest version on your home directory by following the procedure stated above.
• Run each of the following example program but make sure you a re using you installed version. Show to the instructor the compilation line.

Important NOTE: for g++, use the prefix -9 in the configure line to put that as suffix to the commands and avoid collisions with the compiler already installed in the system. This can be done by adding the flag --program-suffix=-9 to the configure command.

Spack was designed by llnl and is targeted to simplify the installation and managment of HPC programs with many possible versions and dependencies. Check the docs at https://spack.readthedocs.io/en/latest/tutorial.html or https://bluewaters.ncsa.illinois.edu/webinars/software-ecosystems/spack . For now let’s just install it, configure it, and install some simple tools.

Fix the following code

#include <iostream>
#include <cstdlib>

int main(int argc, char **argv)
{
  // declare variables
  const int N = 10;
  double * array;
  //reserve memory
  array = new double [N];
  // initialise array
  for (int ii = 0; ii < N; ++ii) {
    array[ii] = 2*(++ii); // != 2*(ii++) ?
  }
  // print some values
  int idx = 2;
  std::cout << array[idx] << std::endl;
  idx = 10;
  std::cout << array[idx] << std::endl; // compiles, but ...

  // free memory
  //delete [] array;

  return EXIT_SUCCESS;
}

Besides the typical printing to help debugging, a first helping hand would the compiler sanitizers. Try recompiling the previous code using the sanitizer flags as

g++ -fsanitize=address -fsanitize=leak -fsanitize=undefined source.cpp

and then run it

./a.out

Fo you observe something new? sanitizers check your code in runtime and help you find errors faster.

The Data Display Debugger - ddd is a graphical frontend for command line debuggers like gdb, dbx, jdb, bashdb, etc, and some modified python debugger. It easies the debug processing by adding some capabilities, like visual monitoring of data and data plotting (killer feature).

NOTE: DDD is hard to configure, it is so rebel. If after configuring some options you are having issues, like dd hanging with a message “Waiting for gdb to start”, or the plotting window is always waiting, then close ddd, delete the config file with rm -rf ~/.ddd, and start over by configuring ONLY one option at a time (configure one option, then close it, then open, then configure other option, then close, etc).

Let’s use ddd to explore some other exercises.

Valgrind official web page

From Web Page: > Valgrind is an instrumentation framework for building dynamic analysis tools. There are Valgrind tools that can automatically detect many memory management and threading bugs, and profile your programs in detail. You can also use Valgrind to build new tools.

The Valgrind distribution currently includes six production-quality tools: a memory error detector, two thread error detectors, a cache and branch-prediction profiler, a call-graph generating cache and branch-prediction profiler, and a heap profiler. It also includes three experimental tools: a heap/stack/global array overrun detector, a second heap profiler that examines how heap blocks are used, and a SimPoint basic block vector generator. It runs on the following platforms: X86/Linux, AMD64/Linux, ARM/Linux, PPC32/Linux, PPC64/Linux, S390X/Linux, MIPS/Linux, ARM/Android (2.3.x and later), X86/Android (4.0 and later), X86/Darwin and AMD64/Darwin (Mac OS X 10.6 and 10.7, with limited support for 10.8).

Modularize one of the previous codes into a header a implementation file and a main file. This will help when writing tests.

Our goal here is to learn to use catch2 to test a very simple function extracted from their tutorial. Later we will modularize the code to practice that and write a useful Makefile.

Google test is a famous and advance unit framework that goes well beyond of what is shown here. You are invited to follow the docs to learn more.

Profiling allows you to learn how the computation time was spent and how the data flows in your program. This can be achieved by just printing the time spent on each section using some internal timers, or, better, by using a tool called “profiler” that shows you how time was spent in your program and which functions called which other functions.

Using a profiler is something that comes very handy when you want to verify that your program does what you want it to do, especially when it is not easy to analyze (it may contain many functions or many calls to different functions and so on).

Remember, WE ARE SCIENTISTS, so we want to profile code to optimize the computation time taken by our program. The key idea then becomes finding where (which function/subroutine) is the computation time spent and attack it using optimization techniques (studied in a previous session of this course, check this).

The first approach is to just add timers to your code. This is a good practice and it is useful for a code to report the time spent on different parts. The following example is extracted from here and modified to make it a bit simpler:

#include <cstdio>
#include <cstdlib>

/*
  Tests cache misses.
*/

int main(int argc, char **argv)
{
  if (argc < 3){
    printf("Usage: cacheTest sizeI sizeJ\nIn first input.\n");
    return 1;
  }
  long sI = atoi(argv[1]);
  long sJ = atoi(argv[2]);

  printf("Operating on matrix of size %ld by %ld\n", sI, sJ);

  long *arr = new long[sI*sJ]; // double array.

  // option 1
  for (long i=0; i < sI; ++i)
    for (long j=0; j < sJ; ++j)
      arr[(i * (sJ)) + j ] = i;

  // option 2
  for (long i=0; i < sI; ++i)
    for (long j=0; j < sJ; ++j)
      arr[(j * (sI)) + i ] = i;

  // option 3
  for (int i=0; i < sI*sJ; ++i) arr[i] = i;

  printf("%ld\n", arr[0]);

  return 0;
}

Now let’s add some timers using the chrono c++11 library (you can see more eamples here and here .

#include <cstdio>
#include <cstdlib>
#include <chrono>
#include <iostream>

/*
   Tests cache misses.
*/

template <typename T>
void print_elapsed(T start, T end );

int main(int argc, char **argv)
{
  if (argc < 3){
    printf("Usage: cacheTest sizeI sizeJ\nIn first input.\n");
    return 1;
  }
  long sI = atoi(argv[1]);
  long sJ = atoi(argv[2]);

  printf("Operating on matrix of size %ld by %ld\n", sI, sJ);

  auto start = std::chrono::steady_clock::now();
  long *arr = new long[sI*sJ]; // double array.
  auto end = std::chrono::steady_clock::now();
  print_elapsed(start, end);

  // option 1
  start = std::chrono::steady_clock::now();
  for (long i=0; i < sI; ++i)
    for (long j=0; j < sJ; ++j)
      arr[(i * (sJ)) + j ] = i;
  end = std::chrono::steady_clock::now();
  print_elapsed(start, end);

  // option 2
  start = std::chrono::steady_clock::now();
  for (long i=0; i < sI; ++i)
      for (long j=0; j < sJ; ++j)
         arr[(j * (sI)) + i ] = i;
  end = std::chrono::steady_clock::now();
  print_elapsed(start, end);

  // option 3
  start = std::chrono::steady_clock::now();
  for (int i=0; i < sI*sJ; ++i) arr[i] = i;
  end = std::chrono::steady_clock::now();
  print_elapsed(start, end);

  printf("%ld\n", arr[0]);

  return 0;
}

template <typename T>
void print_elapsed(T start, T end )
{
  std::cout << "Elapsed time in ms: "
        << std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count()
        << "\n";
}

Analyze the results.

There are many types of profilers from many different sources, commonly, a profiler is associated to a compiler, so that for example, GNU (the community around gcc compiler) has the profiler ’gprof’, intel (corporation behind icc) has iprof, PGI has pgprof, etc. Valgrind is also a useful profiler through the cachegrind tool, which has been shown at the debugging section.

This crash course (more like mini tutorial) will focus on using gprof. Note that gprof supports (to some extend) compiled code by other compilers such as icc and pgcc. At the end we will briefly review perf, and the google performance tools.

FACT 1: According to Thiel, “gprof … revolutionized the performance analysis field and quickly became the tool of choice for developers around the world …, the tool is still actively maintained and remains relevant in the modern world.” (from Wikipedia).

FACT 2: Top 50 most influential papers on PLDI (from Wikipedia).

Inthe following we will use the following code as example: Example 1: Basic gprof (compile with gcc)

//test_gprof.c
#include<stdio.h>

void new_func1(void);
void func1(void);
static void func2(void);
void new_func1(void);

int main(void)
{
  printf("\n Inside main()\n");
  int i = 0;

  for(;i<0xffffff;i++);
  func1();
  func2();

  return 0;
}
void new_func1(void);

void func1(void)
{
  printf("\n Inside func1 \n");
  int i = 0;

  for(;i<0xffffffff;i++);
  new_func1();

  return;
}

static void func2(void)
{
  printf("\n Inside func2 \n");
  int i = 0;

  for(;i<0xffffffaa;i++);
  return;
}

void new_func1(void)
{
  printf("\n Inside new_func1()\n");
  int i = 0;

  for(;i<0xffffffee;i++);

  return;
}

Para usar gprof:

1. Compilar con -pg

     gcc -Wall -pg -g test_gprof.c -o test_gprof
   

2. Ejecutar normalmente. Esto crea o sobreescribe un archivo gmon.out

   ./test_gprof
   

3. Abrir el reporte usando

     gprof test_gprof gmon.out > analysis.txt
   

This produces a file called ’analysis.txt’ which contains the profiling information in a human-readable form. The output of this file should be something like the following:

  Flat profile:

  Each sample counts as 0.01 seconds.
    %   cumulative   self              self     total
   time   seconds   seconds    calls   s/call   s/call  name
   39.64      9.43     9.43        1     9.43    16.79  func1
   30.89     16.79     7.35        1     7.35     7.35  new_func1
   30.46     24.04     7.25        1     7.25     7.25  func2
    0.13     24.07     0.03                             main

   %         the percentage of the total running time of the
  time       program used by this function.

  cumulative a running sum of the number of seconds accounted
   seconds   for by this function and those listed above it.

   self      the number of seconds accounted for by this
  seconds    function alone.  This is the major sort for this
             listing.

  calls      the number of times this function was invoked, if
             this function is profiled, else blank.

   self      the average number of milliseconds spent in this
  ms/call    function per call, if this function is profiled,
         else blank.

   total     the average number of milliseconds spent in this
  ms/call    function and its descendents per call, if this
         function is profiled, else blank.

  name       the name of the function.  This is the minor sort
             for this listing. The index shows the location of
         the function in the gprof listing. If the index is
         in parenthesis it shows where it would appear in
         the gprof listing if it were to be printed.
  
               Call graph (explanation follows)


  granularity: each sample hit covers 2 byte(s) for 0.04% of 24.07 seconds

  index % time    self  children    called     name
                                                   <spontaneous>
  [1]    100.0    0.03   24.04                 main [1]
                  9.43    7.35       1/1           func1 [2]
                  7.25    0.00       1/1           func2 [4]
  -----------------------------------------------
                  9.43    7.35       1/1           main [1]
  [2]     69.7    9.43    7.35       1         func1 [2]
                  7.35    0.00       1/1           new_func1 [3]
  -----------------------------------------------
                  7.35    0.00       1/1           func1 [2]
  [3]     30.5    7.35    0.00       1         new_func1 [3]
  -----------------------------------------------
                  7.25    0.00       1/1           main [1]
  [4]     30.1    7.25    0.00       1         func2 [4]
  -----------------------------------------------

   This table describes the call tree of the program, and was sorted by
   the total amount of time spent in each function and its children.

   Each entry in this table consists of several lines.  The line with the
   index number at the left hand margin lists the current function.
   The lines above it list the functions that called this function,
   and the lines below it list the functions this one called.
   This line lists:
       index  A unique number given to each element of the table.
          Index numbers are sorted numerically.
          The index number is printed next to every function name so
          it is easier to look up where the function is in the table.

       % time This is the percentage of the `total' time that was spent
          in this function and its children.  Note that due to
          different viewpoints, functions excluded by options, etc,
          these numbers will NOT add up to 100%.

       self   This is the total amount of time spent in this function.

       children   This is the total amount of time propagated into this
          function by its children.

       called This is the number of times the function was called.
          If the function called itself recursively, the number
          only includes non-recursive calls, and is followed by
          a `+' and the number of recursive calls.

       name   The name of the current function.  The index number is
          printed after it.  If the function is a member of a
          cycle, the cycle number is printed between the
          function's name and the index number.


   For the function's parents, the fields have the following meanings:

       self   This is the amount of time that was propagated directly
          from the function into this parent.

       children   This is the amount of time that was propagated from
          the function's children into this parent.

       called This is the number of times this parent called the
          function `/' the total number of times the function
          was called.  Recursive calls to the function are not
          included in the number after the `/'.

       name   This is the name of the parent.  The parent's index
          number is printed after it.  If the parent is a
          member of a cycle, the cycle number is printed between
          the name and the index number.

   If the parents of the function cannot be determined, the word
   `<spontaneous>' is printed in the `name' field, and all the other
   fields are blank.

   For the function's children, the fields have the following meanings:

       self   This is the amount of time that was propagated directly
          from the child into the function.

       children   This is the amount of time that was propagated from the
          child's children to the function.

       called This is the number of times the function called
          this child `/' the total number of times the child
          was called.  Recursive calls by the child are not
          listed in the number after the `/'.

       name   This is the name of the child.  The child's index
          number is printed after it.  If the child is a
          member of a cycle, the cycle number is printed
          between the name and the index number.

   If there are any cycles (circles) in the call graph, there is an
   entry for the cycle-as-a-whole.  This entry shows who called the
   cycle (as parents) and the members of the cycle (as children.)
   The `+' recursive calls entry shows the number of function calls that
   were internal to the cycle, and the calls entry for each member shows,
   for that member, how many times it was called from other members of
   the cycle.

  
  Index by function name

     [2] func1                   [1] main
     [4] func2                   [3] new_func1

The output has two sections: >*Flat profile:* The flat profile shows the total amount of time your program spent executing each function.

Call graph: The call graph shows how much time was spent in each function and its children.

See https://www.brendangregg.com/FlameGraphs/cpuflamegraphs.html

Valgrind allows not only to debug a code but also to profile it. Here we will see how to use cachegrind, to check for cache misses, and callgrind, for a calling graph much like tools like perf and gprof.

See https://github.com/gperftools/gperftools

From Google performance tools : These tools are for use by developers so that they can create more robust applications. Especially of use to those developing multi-threaded applications in C++ with templates. Includes TCMalloc, heap-checker, heap-profiler and cpu-profiler.

In brief:

  TC Malloc:

  gcc [...] -ltcmalloc
  Heap Checker:

  gcc [...] -o myprogram -ltcmalloc
  HEAPCHECK=normal ./myprogram
  Heap Profiler:

  gcc [...] -o myprogram -ltcmalloc
  HEAPPROFILE=/tmp/netheap ./myprogram
  Cpu Profiler:

  gcc [...] -o myprogram -lprofiler
  CPUPROFILE=/tmp/profile ./myprogram

Basically, when ypu compile, you link with the required library. Then, you can generate a callgraph with profiler info. Try to install the google performance tool on your home and test them with the previous codes. Please review the detailed info for each tool: for example, for the cpu profiler, check Cpu profiler info

1. Take this code and modify it (separate its behaviour into functions) to be able to profile it.
2. Download https://bitbucket.org/iluvatar/scientific-computing-part-01/downloads/CodigosIvan.zip, run , and optimize them.
3. Experiment: Take this code, profile it and try to optimize it in any way.

See : https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html

Sometimes just using compiler optimization flags with can improve dramatically our code performance. Normally, you want to run your code through several compiler flags to check for speed ups but also for possible uncovered bugs.

Compile without (-O0) and with optimization (-O3) the following code and compare their runtime:

// Credits : Ivan Pulido
/* Shows a way to do operations that require a specific order (e.g.,
 * transpositions) while avoiding cache misses. */

#include <stdio.h>
#include <time.h>
#include <stdlib.h>

#define min( a, b ) ( ((a) < (b)) ? (a) : (b) )

int main(){
    const int n = 512;
    const int csize = 32;
    float ***a, ***b;
    clock_t cputime1, cputime2;
    int i,j,k,ii,jj,kk;

    // Allocating memory for array/matrix
    a = malloc(n*sizeof(float **));
    for (i=0; i<n; i++){
        a[i] = malloc(n*sizeof(float*));
        for (j=0; j<n; j++)
            a[i][j] = malloc(n*sizeof(float));
    }
    b = malloc(n*sizeof(float **));
    for (i=0; i<n; i++){
        b[i] = malloc(n*sizeof(float*));
        for (j=0; j<n; j++)
            b[i][j] = malloc(n*sizeof(float));
    }

    // Filling matrices with zeros
    for(i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            for (k=0; k<n; ++k)
                a[i][j][k] = 0;
    for(i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            for (k=0; k<n; ++k)
                b[i][j][k] = 0;

    // Direct (inefficient) transposition
    cputime1 = clock();
    for (i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            for (k=0; k<n; ++k)
                a[i][j][k] = b[k][j][i];
    cputime2 = clock() - cputime1;
    printf("Time for transposition: %f\n", ((double)cputime2)/CLOCKS_PER_SEC);

    // Transposition using cache-blocking
    cputime1 = clock();
    for (ii=0; ii<n; ii+=csize)
        for (jj=0; jj<n; jj+=csize)
            for (kk=0; kk<n; kk+=csize)
                for (i=ii; i<min(n,ii+csize-1); ++i)
                    for (j=jj; j<min(n,jj+csize-1); ++j)
                        for (k=kk; k<min(n,kk+csize-1); ++k)
                            a[i][j][k] = b[k][j][i];
    cputime2 = clock() - cputime1;
    printf("Time for transposition: %f\n", ((double)cputime2)/CLOCKS_PER_SEC);

    return 0;
}

When you compile without optimization and execute, you will get something like

gcc -O0 cache_blocking.c
./a.out

But, if you use optimization, you can get an important improvement

gcc -O2 cache_blocking.c
./a.out

Actually we can compare all optimization levels to check their impact:

for level in 0 1 2 3; do echo "level: $level"; gcc -O$level cache_blocking.c; ./a.out; done

Be careful when using -O3: it activates some unsafe math optimizations that will not respect IEEE754.

The memory layout should always be taken into account since cache misses will greatly affect a program performance. Always measure with a profiler and if possible with cachegrind in order to detect possible excessive cache misses.

The following code shows how a simple index change in a matrix operation could have a huge impact on the app performance:

// Credit Ivan Pulido
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(){
    const int n = 256;
    clock_t cputime1, cputime2;
    float ***a;
    int i,j,k;

    // Allocating memory for array/matrix
    a = malloc(n*sizeof(float **));
    for (i=0; i<n; i++){
        a[i] = malloc(n*sizeof(float*));
        for (j=0; j<n; j++)
            a[i][j] = malloc(n*sizeof(float));
    }
    cputime1 = clock();
    for (k=0; k<n; ++k)
        for (j=0; j<n; ++j)
            for (i=0; i<n; ++i)
                a[i][j][k] = 1.0;
    cputime2=clock() - cputime1;
    printf("Time with fast index inside: %lf\n", ((double)cputime2)/CLOCKS_PER_SEC);

    cputime1 = clock();
    for(i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            for (k=0; k<n; ++k)
                a[i][j][k] = 2.3;
    cputime2=clock() - cputime1;
    printf("Time with fast index outside: %lf\n", ((double)cputime2)/CLOCKS_PER_SEC);

    // Clearing memory
    for (i=0; i<n; i++){
        for (j=0; j<n; j++)
            free(a[i][j]);
        free(a[i]);
    }
    free(a);

    return 0;
}

gcc codes/cache_lines.c
./a.out

Being nice to the cache comes from the modern cpu architecture, as shown in the next figure

The relative bandwidth across the different cpu and momory controllers explain why it is important to have the processed data as much time as possible in the cache

Since c/c++ store matrices in row major order, then accesing the memory in the same way will benefit the cache and therefore increases our program performance.

Scientific libraries are written by people with deep knowledge of the computer and algorithms internals , therefore saving developer time and resources. One should always try to use a established scientific library instead of writting everything from scratch. Sometimes even automatic parallelization comes for free. Examples are the Intel or AMD math libraries, the Cuda toolkit for programming on nvidia cards, and so on.

The following examples shows the time taken to transpose a matrix using a traditional approach, a blocking approach , and the eigen c++ library. giving the following results:

Time for memory allocation: 0.315618
Time for filling: 3.607817
Time for direct transposition: 2.870691
Time for blocked transposition: 0.380954
Time for transposition with eigen: 0.000003
Time for transposition with copy in eigen: 0.000001
Time for transposition with full copy in eigen: 3.339495
-2999.9
-2999.9
15200.1
15200.1
15200.1

This shows, for instance, that eigen in some cases is not even computing the transpose but just creating an expression to access the original matrix, hence the huge speeedup and also the slow down when fully creating a copy.

There are other techniques sometimes applied in very speficic situations, like

• Loop interchanges
• Loop unrolling

for(int ii = 0; ii < n; ii++) {
    array[ii] = 2*ii;
}

for(int ii = 0; ii < n; ii += 3) {
    array[ii] = 2*ii;
    array[ii+1] = 2*(ii+1);
    array[ii+2] = 2*(ii+2);
}

• Loop Fusion/Fision
• Prefetching
• Floating point division
• Vectorization
• etc

In general those techniques can applied after a careful determination that they are really needed, and sometimes the compilers already apply them at some optimization levels. Therefore it is advisable to focus on the clarity of the program first and let the compiler do its job.

• https://www.agner.org/optimize/

https://www.gnu.org/software/gsl/ Documentation: https://www.gnu.org/software/gsl/doc/html/index.html

Check the standard library for implementatios of many special functions. Must compile with -std=c++17

http://eigen.tuxfamily.org/index.php?title=Main_Page Start by: http://eigen.tuxfamily.org/dox/GettingStarted.html

• Basic matrix

 1: #include <iostream>
 2: #include <eigen3/Eigen/Dense>
 3: 
 4: int main(void)
 5: {
 6:     Eigen::MatrixXd m(2, 2);
 7:     m(0, 1) = -9.877;
 8:     std::cout << m << std::endl;
 9: 
10:     return 0;
11: }
12: 

• Exercise: Create a 8x8 matrix and print it to the screen.
• random matrix

#include <iostream>
#include <cstdlib>
#include <eigen3/Eigen/Dense>

int main(int argc, char **argv)
{
    //const int SEED = std::atoi(argv[1]);
    //srand(SEED); // control seed
    Eigen::MatrixXd m = Eigen::MatrixXd::Random(5,5);
    std::cout << m << std::endl;

    return 0;
}

• Exercise: Create a 8x8 random matrix and print it to the screen. Control the seed, print it twice to check that the elements are the same.
• First example from the tutorial:

  // First example on how to use eigen
  // compile with: g++ -std=c++11 -I /usr/include/eigen3
  #include <iostream>
  #include <eigen3/Eigen/Dense>

  int main()
  {
    // X = Dynamic size
    // d = double
    Eigen::MatrixXd m(2,2);
    m(0,0) = 3;
    m(1,0) = 2.5;
    m(0,1) = -1;
    m(1,1) = m(1,0) + m(0,1);
    std::cout << m << std::endl;
  }

• Second example from the tutorial: Matrix-vector multiplication

#include <iostream>
#include <eigen3/Eigen/Dense>

int main()
{
    const int N = 3;
    srand(10);
    Eigen::MatrixXd m = Eigen::MatrixXd::Random(N, N);
    m = (m + Eigen::MatrixXd::Constant(N, N, 1.2)) * 50;
    std::cout << "m =" << std::endl << m << std::endl;
    Eigen::VectorXd v(N);
    v << 1, 2, 3;
    std::cout << "m * v =" << std::endl << m * v << std::endl;
}

In the following we want to explore the performance of the following program as a function of the system size

#include <iostream>
#include <eigen3/Eigen/Dense>

void solvesystem(int size);

int main(int argc, char ** argv)
{
  const int M = atoi(argv[1]); // Matrix size
  solvesystem(M);
}

void solvesystem(int size)
{
  Eigen::MatrixXd A = Eigen::MatrixXd::Random(size, size);
  Eigen::MatrixXd b = Eigen::MatrixXd::Random(size, 1);
  // crono star
  //Eigen::MatrixXd x = A.fullPivLu().solve(b);
  Eigen::MatrixXd x = A.lu().solve(b);
  // crono end
  double relative_error = (A*x - b).norm() / b.norm(); // norm() is L2 norm
  std::cout << A << std::endl;
  std::cout << b << std::endl;
  std::cout << "The relative error is:\n" << relative_error << std::endl;
}

// usar chrono para medir el tiempo
// gprof
// remove printing
// modify printing

#include <iostream>
#include <eigen3/Eigen/Dense>

void solvesystem(int size);

int main(int argc, char ** argv)
{
  const int M = atoi(argv[1]); // Matrix size
  solvesystem(M);
}

void solvesystem(int size)
{
  Eigen::MatrixXd A = Eigen::MatrixXd::Random(size, size);
  Eigen::MatrixXd b = Eigen::MatrixXd::Random(size, 1);
  // crono star
  Eigen::MatrixXd x = A.fullPivLu().solve(b);
  // crono end
  // double relative_error = (A*x - b).norm() / b.norm(); // norm() is L2 norm
  //std::cout << A << std::endl;
  //std::cout << b << std::endl;
  //std::cout << "The relative error is:\n" << relative_error << std::endl;
}

// usar chrono para medir el tiempo
// gprof
// remove printing
// modify testing

If you use just

time ./a.out

you will measuring the whole execution time, including the memory allocation, the actual operation, and the memory freeing. It is much better to actually measure only what we are interesting in, which is the actual memory operation.

To create a reasonable benchmark, we need to

• Fill the matrix randomly, controlling the seed
• Measure only processing time (wall and cpu times)
• Statistics
• Print the time as function of N and plot (run this in the shell)
• Use threads (openmp)

Now do the same study but with matrix multiplication. Implement the matrix matmul with eigen (see the code snippet at the end). Always compile with the flags “-fopenmp -O3”. And now execute in the following way

OMP_NUM_THREADS=VAR ./a.out ...

where VAR takes values from 1, 2, 3, 4, … up the output from the command nproc. This means that you are starting to use parallelization. You must run this in a multicore system (or give your virtual machine enough resources) with a least 8 threads. Plot the time (with error bars) as a function of the system size, and you will have four curves on that graph, each one for each VAR value.

Basic snippet for matmul in eigen:

double multiply(int size)
{
  // create matrices
  Eigen::MatrixXd A = Eigen::MatrixXd::Random(size, size);
  Eigen::MatrixXd B = Eigen::MatrixXd::Random(size, size);

  auto start = std::chrono::steady_clock::now();
  auto C{A*B}; // MULTIPLY
  double tmp = C(0, 0); // use the matrix to make eigen compute it
  auto end = std::chrono::steady_clock::now();

  std::clog << tmp << std::endl; // use the matrix to make eigen compute it

  std::chrono::duration<double> elapsed_seconds = end-start;
  return elapsed_secons.count();
}

See presentations hpc-intro.pdf, 07-HPC-Distributed.pdf and 07-Parallel-general-metrics.pdf . Special care to metrics.

• https://hpc.llnl.gov/documentation/tutorials

Debuggers:

• http://www.cs.uoregon.edu/research/tau/home.php
• https://vampir.eu/

There is a point where a serial version of our code is not the most optimal way to exploit our computational resources (but it might be in the case of embarrassingly parallel problems where you can just run several programs at once). For instance, you might want to use all the cores on your multicore system, ideally reducing the execution time, or you need to explore larger system sizes that could consume a lot of memory or need too much time.

Typically, Moore’s law allowed to wait for a bit in order to get a better machine so your algorithms will run faster.

But due to physics limitation and power considerations, it is now typical to have multicore systems

Recently, power considerations are being more and more relevant:

• https://en.wikipedia.org/wiki/Performance_per_watt?useskin=vector
• https://en.wikipedia.org/wiki/Koomey%27s_law?useskin=vector
• https://www.apple.com/newsroom/2022/03/apple-unveils-m1-ultra-the-worlds-most-powerful-chip-for-a-personal-computer/

• https://www.extremetech.com/extreme/328541-the-apple-m1-pro-and-m1-maxs-power-efficiency-should-rattle-intel-amd)

  

  Figure 6: Power comparison between M1 Max and Intel i9

At the same time, the computational problems size and/or complexity has been steadily increasing in time, requiring distributed computing techniques.

(see also https://en.wikipedia.org/wiki/Flynn%27s_taxonomy?useskin=vector). Recently, besides CPU parallelization, the GPU parallelization has become very relevant (see https://en.wikipedia.org/wiki/General-purpose_computing_on_graphics_processing_units?useskin=vector ), where CUDA , OpenACC, and others, are the relevant technologies.

In our case, we will be more focused on the cluster computing aspect, while there are more HPC approaches, like grid computing, cloud computing, and so on. One of the goals of HPC is to get better results faster and/or to exploit better current or future resources.

But, as usual, you should always measure. All programs have a serial part that cannot be parallelized and a parallel part than can. Using more processors/threads can reduce only the parallel, so a 100% serial program cannot really take advantage of a parallel system. This is known as Amdahls law, https://en.wikipedia.org/wiki/Amdahl%27s_law?useskin=vector

At the end, the user must also gauge its application performance. Blindly reserve of HPC resources represent a non efficient cluster use, and higher costs. In this regard, parallel metrics are really crucial. The next to figures show the speedup and the parallel efficiency. As you can see, they are limited by the hardware (and algorithms)

• Speedup:
• Parallel efficiency:

There are many aspects to take into account in the HPC field. If you are a user, you should know abount the type of parallelization (shared memory, disitributed memory, gpu programming), the type of hardware you are using, the resource manager, the data storage and so on. The goal of a system administrator is to make that easier, but that is not always possible. Check the 12-12-hkhlr_quick_reference-goethe-hlr (from https://csc.uni-frankfurt.de/wiki/doku.php?id=public:start) for an example of a typical cluster config and offerings.

These are examples from the Archer cluster at https://www.archer2.ac.uk/