% SIAM Article Template
\documentclass[]{article}

\usepackage{arxiv}

\usepackage{xcolor}
\definecolor{darkgreen}{rgb}{0.0, 0.5, 0.0}
\definecolor{amaranth}{rgb}{0.9, 0.17, 0.31}
\definecolor{azure}{rgb}{0.0, 0.5, 1.0}

\usepackage{listings}

\lstset{language=C++,
 basicstyle=\ttfamily,
 keywordstyle=\color{azure}\ttfamily,
 stringstyle=\color{amaranth}\ttfamily,
 commentstyle=\color{darkgreen}\ttfamily,
 morecomment=[l][\color{magenta}]{\#},
 breaklines=true
}

\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{xfrac}

\usepackage{hyperref}

% Information that is shared between the article and the supplement
% (title and author information, macros, packages, etc.) goes into
% ex_shared.tex. If there is no supplement, this file can be included
% directly.


% Optional PDF information
%\ifpdf
%\hypersetup{
%  pdftitle={Advanced Parallel Programming in C++},
%  pdfauthor={P. Diehl}
%}
%\fi

% The next statement enables references to information in the
% supplement. See the xr-hyperref package for details.

%% Use \myexternaldocument on Overleaf


% FundRef data to be entered by SIAM
%<funding-group>
%<award-group>
%<funding-source>
%<named-content content-type="funder-name"> 
%</named-content> 
%<named-content content-type="funder-identifier"> 
%</named-content>
%</funding-source>
%<award-id> </award-id>
%</award-group>
%</funding-group>
\newcommand{\todo}[1]{{~\textbf{\color{red}[TODO: #1]}}}

\begin{document}

\title{Advanced Parallel Programming in C++}

\author{Patrick Diehl, Steven R. Brandt, and Hartmut Kaiser \\
Center of Computation \& Technology \\
Louisiana State University}

\maketitle

% REQUIRED
\begin{abstract}
 Many applied mathematics codes are based on the C++ programming language. However, many of these codes do not utilize modern C++ features, especially the features from the recent C++ 17 standard for parallel computations. These features can significantly simplify parallel computations using multiple cores. This review briefly introduces asynchronous programming facilities introduced in the C++ 11 standard and the parallel algorithms introduced in the C++ 17 standard. With these two paradigms, C++ code can be accelerated using multiple cores using the C++ standard without any need for some external special purpose libraries, like \emph{e.g.} OpenMP. The C++20 standard made additional and valuable contributions to the C++ api for parallelism and concurrency. While these may still be unavailable in some compilers, they are already implemented in HPX~\cite{kaiser2020hpx}.
\end{abstract}

% REQUIRED
%\begin{keywords}
%  example, \LaTeX
%\end{keywords}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Asynchronous programming}
\label{sec:async}
The C++11 standard~\cite{cxx11_standard} has introduced asynchronous task programming. Here, the concept of futurization using \lstinline{std::future} was introduced. A future is a placeholder for the result of an asynchronous function (launched e.g. using \lstinline{std::async()}). The runtime promises that once the computation of the task passed to \lstinline{std::async()} has finished, the future will receive the result. One can use the function \lstinline{.get()} to access the result. If the computation has finished the result will be returned. If the computation is still ongoing, the function will suspend the current thread until the future has become ready, at which point the result is returned to the user. Let us look at the following Taylor expansion for the natural logarithm
\begin{align}
    \ln(x) = \sum\limits_{n=1}^N \frac{(-1)^{n+1}}{n}(x-1)^n \textbf{.}
\end{align}
Listing~\ref{lst:async:cpp} shows the asynchronous computation of the Taylor series. First, a function \lstinline{taylor_part} that computes a part of the sum is defined in Line~\ref{lst:async:cpp:1}. In Line~\ref{lst:async:cpp:2} the sum of length $n$ is separated into two partitions of $\sfrac{n}{2}$ and each of them is asynchronously launched, and a future is returned. Since these two function launches run concurrently on two cores, we need to synchronize their execution and collect the results. In Line~\ref{lst:async:cpp:3} the \lstinline{.get()} function is called to access the result of each of the futures. Note that the code will introduce a barrier and waits until the first and second future are ready.


\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Asynchronous computation of the Taylor expansion for the natural logarithm.,
 label={lst:async:cpp},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <cmath>
#include <future>
#include <iostream>

double taylor_part(double x, size_t start, size_t end){ |\label{lst:async:cpp:1}|
    
    double result = 0; 
    for(size_t i = start; i < end; i++)
    {
        result += std::pow(-1,i+1)/i * std::pow(x-1,i);
    }
    return result;
}

int main(void){

    double x = 0.5;
    size_t n = 100;

    // Launch the function async |\label{lst:async:cpp:2}|
    std::future<double> f1 = std::async(taylor_part,x,1,n/2);
    std::future<double> f2 = std::async(taylor_part,x,n/2,n);

    // Gather the result 
    std::cout << f1.get() + f2.get() << std::endl;  |\label{lst:async:cpp:3}|
    return EXIT_SUCCESS;
}
\end{lstlisting}

The C++ standard library for parallelism and concurrency (HPX)~\cite{kaiser2020hpx} implements the same functionality as mandated by the C++ standard. However, additional, more advanced synchronization facilities are provided by HPX. The function \lstinline{.then()} provides a function similar to \lstinline{.get()}, however it takes a function that is only to be run after the future's value arrives. While slightly less intuitive to program, this method avoids blocking. Also, it is possible to use \lstinline{when_all()}, which waits for one or more futures. Once all futures are ready, \lstinline{when_all()} returns a future that becomes ready only after all argument futures have become ready. Listing~\ref{lst:async:hpx} uses these features.


\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Asynchronous computation of the Taylor expansion for the natural logarithm using HPX.,
 label={lst:async:hpx},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <hpx/hpx_main.hpp>   
#include <hpx/future.hpp>
#include <cmath>
#include <iostream>

double taylor_part(double x, size_t start, size_t end){
    
    double result = 0; 
    for(size_t i = start; i < end; i++)
    {
        result += std::pow(-1,i+1)/i * std::pow(x-1,i);
    }
    return result;
}

int main(void){

    double x = 0.5;
    size_t n = 100;
    size_t partitions = 3;
    std::vector<hpx::future<double>> futures;
    futures.reserve(partitions);
    
    // Launch the function async
    futures.push_back(hpx::async(taylor_part,x,1,n/3));
    for(size_t i = 1; i < partitions; i++)
      futures.push_back(hpx::async(taylor_part,x,i*n/3,(i+1)*n/3));

    // Gather the result
    hpx::when_all(futures).then([](auto&& f){
        std::vector<hpx::future<double>> futures = f.get();
        double result = 0;
        for(size_t i = 0; i < futures.size(); i++)
            result += futures[i].get();
        std::cout << result << std::endl;
        });
    return EXIT_SUCCESS;
}
\end{lstlisting}

Since the C++20 standard for coroutines has emerged, HPX futures have become awaitable. That means that the less natural style of programming using \lstinline{.then()} can be replaced by calls to \lstinline{co_await}. Listing~\ref{lst:coawait:hpx} illustrates how to use this style of programming with HPX. The call to \lstinline{co_await} suspends execution of the current thread until the value in the future is ready. Therefore, \lstinline{co_await} can return the value directly.

\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Asynchronous computation of the Taylor expansion for the natural logarithm using HPX.,
 label={lst:coawait:hpx},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <hpx/hpx_main.hpp>   
#include <hpx/future.hpp>
#include <cmath>
#include <iostream>

double taylor_part(double x, size_t start, size_t end){
    
    double result = 0; 
    for(size_t i = start; i < end; i++)
    {
        result += std::pow(-1,i+1)/i * std::pow(x-1,i);
    }
    return result;
}

// Can only use co_await inside a function that returns a future
hpx::future<double> get_taylor_part(double x, size_t n, size_t partitions) {
    double result = 0;
    std::vector<hpx::future<double>> futures;
    futures.reserve(partitions);
    
    // Launch the function async
    futures.push_back(hpx::async(taylor_part,x,1,n/3));
    for(size_t i = 1; i < partitions; i++)
      futures.push_back(hpx::async(taylor_part,x,i*n/3,(i+1)*n/3));

    // Gather the results
    auto futures2 = co_await hpx::when_all(futures);

    // Sum the results
    for(size_t i=0; i < futures2.size(); i++)
        result += co_await futures2[i];

    co_return result;
}

int main(void){
    double x = 0.5;
    size_t n = 100;
    size_t partitions = 3;
    std::cout << get_taylor_part(x,n,partitions).get() << std::endl;
    return EXIT_SUCCESS;
}
\end{lstlisting}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parallel algorithms}
\label{sec:parallel}
A common opportunity for shared-memory parallel programming in C++ is Open Multi-Processing (OpenMP)~\cite{dagum1998openmp}. Let us look at the implementation of computing the square root of all elements in a vector $v$ of size $n$ 
\begin{align}
    v_i = \sqrt{v_i}, \; \forall i \in {1,\ldots,n}
\end{align}
in a parallel fashion using OpenMP (see Listing~\ref{lst:openmp}). In Line~\ref{lst:openmp:1}, a \lstinline|std::vector| is generated with length $n=10000$ and in Line~\ref{lst:openmp:2} the vector is filled with random numbers using the algorithms provided by the C++ standard library. In Line~\ref{lst:openmp:5} a \lstinline|for| loop is used to iterate over the elements of the vector. Now, the \lstinline|#pragma|-based OpenMP API is used to parallelize the \lstinline|for| loop. In Line~\ref{lst:openmp:4} the \lstinline|#pragma| to execute the loop in parallel is added. However, an additional API is needed to achieve the shared-memory parallelism. 

\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Example to compute the element-wise square root of a vector using OpenMP.,
 label={lst:openmp},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <cmath>
#include <vector>
#include <numeric>
#include <random>
#include <algorithm>
   
int main(void){
   
    size_t n = 10000;
    
    // Generate the vector with length n
    std::vector<double> v(n); |\label{lst:openmp:1}|
    
    // Initialize the vector with -1
    std::generate(v.begin(),v.end(),std::rand);|\label{lst:openmp:2}|
  
    #pragma omp parallel for |\label{lst:openmp:4}|
    for (size_t i = 0 ; i < v.size(); i++) |\label{lst:openmp:5}|
        v[i] = std::sqrt(i);

    return EXIT_SUCCESS;
}
\end{lstlisting}



In the C++ 17 standard~\cite{cxx17_standard}, parallel algorithms were introduced. The shared-memory parallelism they provide is based on the Thread Building Blocks library (TBB)~\cite{10.5555/1352079.1352134}. While it is possible to write parallel code directly using the C++ standard, and there is no need to use a second API, the C++ standard library contains $69$ algorithms which work on standard containers, e.g. \lstinline|std::vector| or \lstinline|std::list|. Some versions of C++ libraries may not implement all of them. Also in C++ 17, execution policies were introduced to these $69$ algorithms, allowing them to execute in parallel. Depending on the version of your C++ libraries, this feature may also be incomplete.

In this next listing, we will revise the previous example to use parallel algorithms. See Listing~\ref{lst:for:each}. Up to Line~\ref{lst:for:each:1}, the code is unchanged, except for the new header file \lstinline|#include <execution>|, which is needed for the execution policies. In Line~\ref{lst:for:each:2}, a vector containing the indices of each vector element is generated. In Line~\ref{lst:for:each:3} the OpenMP parallel \lstinline|for| loop is replaced by \lstinline|std::for_each| from the C++ standard library algorithms. Note that this was possible before. However, the first argument of this function, the execution policy, defines that this algorithm is to be executed in parallel. The C++ 17 standard introduced the following execution policies:
\begin{itemize}
    \item \textbf{Sequential execution}:\\
    By adding \lstinline|std::execution::seq| the algorithm is executed sequentially using one thread as in the previous C++ standards.
    \item \textbf{Parallel execution}:\\
    By adding \lstinline|std::execution::par| the algorithm is executed in parallel using all available threads. 
    \item \textbf{Vectorized parallel execution}:\\
    By adding \lstinline|std::execution::par_unseq| the algorithm is executed in parallel but in addition vectorization is used.
\end{itemize}
By specifying the execution policy in Line~\ref{lst:for:each:4}, the algorithm is easily parallelized. Note that this was implemented using only the C++ standard and not an external tool like OpenMP.


\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Example: Compute the element-wise square root of a vector using C++.,
 label={lst:for:each},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <cmath>
#include <vector>
#include <numeric>
#include <random>
#include <algorithm>
#include <execution>
   
int main(void){
   
    size_t n = 10000;
    // Generate the vector with length n
    std::vector<double> v(n); 
    // Initialize the vector with -1
    std::generate(v.begin(),v.end(),std::rand);|\label{lst:for:each:1}|
  
    // Fill a new vector containing the indices |\label{lst:for:each:2}|     
    std::vector<size_t> i = std::vector<size_t>(n);
    std::iota(std::begin(i), std::end(i), 0);
    
    std::for_each(|\label{lst:for:each:3}| 
        std::execution::par,|\label{lst:for:each:4}| 
        i.begin(),
        i.end(),
        [&](auto&& item)
        {
            v[item] = std::sqrt(v[item]);
        });

    return EXIT_SUCCESS;
}
\end{lstlisting}

Next, we will look at an example that computes the sum $s$ of all elements in a vector $v$ with $n$ elements
\begin{align}
    s = \sum\limits_{i=0}^n v_i \text{,} 
\end{align}
see Listing~\ref{lst:reduce}. In Line~\ref{lst:reduce} a \lstinline|std::vector| with length \lstinline|n| is generated. The first algorithm of the C++ Standard Library \lstinline|std::fill| is used to fill all elements with negative one, see Line~\ref{lst:reduce:2}. Finally, the second algorithm \lstinline|std::accumulate| is used to compute the sum of all elements. Note that the naive approach would be to use a \lstinline|for| loop to iterate over all the elements. In Line~\ref{lst:reduce:4} the result is printed for validation. Again, one can use the execution policies to easily parrallelize this code by just adding the execution policy in Line~\ref{lst:reduce:par:2} and Line~\ref{lst:reduce:par:3}, respectively. The corresponding code is shown in Listing~\ref{lst:reduce:par}.

\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Example to compute the sum of all elements in the vector using the C++ standard library.,
 label={lst:reduce},
 float=tb,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <iostream>
#include <vector>
#include <numeric>

int main(void){

    size_t n = 10000;
    // Generate the vector with length n
    std::vector<int> v(n); |\label{lst:reduce:1}|
    // Initialize the vector with -1
    std::fill(n2.begin(),n2.end(),-1); |\label{lst:reduce:2}|

    // Compute the sum of all elements 
    int s = std::accumulate(v.begin(),v.end(),0.0); |\label{lst:reduce:3}|

    // Output the result
    std::cout << "Result= " << result << std::endl;|\label{lst:reduce:4}|

    return EXIT_SUCCESS;
}
\end{lstlisting}

\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Example to compute the sum of all elements in the vector using the parallel algorithms provided by the C++ 17 standard.,
 label={lst:reduce:par},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip,
 escapechar=|
 ]
#include <iostream>
#include <vector>
#include <numeric> |\label{lst:reduce:par:1}|
#include <execution>

int main(void){

  constexpr size_t n = 10000;

  // Generate the vector with length n
  std::vector<int> v(n); 

  // Initialize the vector with -1
  std::fill(std::execution::par,v.begin(),v.end(),-1); |\label{lst:reduce:par:2}|

  // Compute the sum of all elements 
  int s = std::reduce(std::execution::par,v.begin(),v.end(),0.0); |\label{lst:reduce:par:3}|

  // Output the result
  std::cout << "Result= " << result << std::endl;

  return EXIT_SUCCESS;
}
\end{lstlisting}

These two small examples show only a few of the parallel algorithms of the C++ 17 standard, which are experimentally supported by the GNU compiler collection (gcc) $\geq$ 9 and the Microsoft Visual Studio compiler (MVSC). However, the C++ standard library for parallelism and concurrency (HPX) implements all the defined features in the C++ 17 standard. Note that you can replace the \lstinline|std::| name space by \lstinline|hpx::| name space, since HPX strictly enhances the C++ 17 standard. We will show two small examples with new features which are not yet implemented in the major implementations. First, HPX provides some easier way to iterate over a vector without generating a vector of indices using \lstinline|std::for_each|, see Listing~\ref{lst:for:each:hpx}. Listing~\ref{lst:for:each:hpx} shows the simplification for iterating over the vector. In the previous example, we had to generate the vector of indices, which increases the memory consumption of the program. In Line~\ref{lst:for:each:hpx:2}, we use a \lstinline|hpx::for_loop| which is more like the \lstinline|for| loop where the start index is provided as the second argument and the end index as the third argument. In the lambda function, we can access the index, $i$ which will go from zero up to \lstinline|v.size()|. In addition, the parallel algorithms can be combined with the concept of futurization to obtain an \lstinline{hpx::future}. Note that this is not provided by the C++ 17 standard.

\begin{lstlisting}[
 basicstyle=\small\ttfamily,
 caption=Example to compute the sum of all elements in the vector using the parallel algorithms provided by HPX.,
 label={lst:for:each:hpx},
 float=htbp,
 numbers=left,
 escapechar=|,
 tabsize=2,
 numberstyle=\tiny,
 numbersep=5pt,
 belowskip=-0.8 \baselineskip
 ]
#include <hpx/hpx_main.hpp>
#include <hpx/algorithm.hpp>
#include <hpx/future.hpp>

#include <iostream>
#include <vector>
#include <numeric>

int main(void) {
  constexpr size_t n = 10000;
  
  // Generate the vector with length n
  std::vector<int> v(n);
  
  // Initialize the vector with -1
  hpx::ranges::fill(hpx::execution::par, v, -1); |\label{lst:for:each:hpx:1}|
  
  // Compute the sum of all elements |\label{lst:for:each:hpx:2}|
  hpx::future<double> f = hpx::ranges::reduce(
      hpx::execution::par(hpx::execution::task), v, 0.0);
  
  // Output the result
  std::cout << "Result= " << f.get() << std::endl;

  return EXIT_SUCCESS;
}
\end{lstlisting}


\section{Conclusions}
\label{sec:conclusions}

In this article we have shown that it is possible to write parallel C++ programs that take advantage of features in the C++17 and C++20 standards, both to provide task parallelism and loop parallelism (with the parallel algorithms). In addition, we've shown how to use more advanced parallel features only available in the HPX programming framework.

Moreover, we have identified two science applications which base their parallelism on this functionality. First, PeriHPX~\cite{jha2021nlmech,diehl2020asynchronous} a non-local continuums mechanics code using all these features for its shared memory implementation. Second, Octo-Tiger~\cite{10.1093/mnras/stab937} a multi-physics astrophysics code for simulating stellar mergers.  

\appendix

\section*{Supplementary materials}
The source code of all the examples is available on GitHub\footnote{\url{https://github.com/diehlpk/modern-cpp-examples}} or Zenodo~\cite{}, respectively. The build system CMake is used to compile all the examples.  

\section*{Acknowledgments}
I would like to acknowledge the students of the course ``Parallel Computational Math'' offered at the Mathematics Department at Louisiana State University for their feedback while teaching these modern C++ features as part of the course.

\bibliographystyle{plain}
\bibliography{references}
\end{document}
