C++ in HPC & Scientific Computing

Why C++ dominates HPC

Scientific computing demands maximum floating-point throughput, minimal overhead, and control over memory layout and hardware features. C++ delivers this while providing abstractions that make complex algorithms expressible. Languages like Python delegate to C++ libraries (NumPy calls BLAS, PyTorch calls CUDA C++).

Eigen — linear algebra

Eigen is the de facto C++ linear algebra library. It uses expression templates to eliminate temporaries and generate code equivalent to hand-tuned BLAS routines — without you writing a single intrinsic.

cpp

Godbolt

#include <Eigen/Dense>

using Matrix4d = Eigen::Matrix4d;
using Vector3d = Eigen::Vector3d;
using MatrixXd = Eigen::MatrixXd;

// Fixed-size — stack-allocated, fully inlined, no heap
Matrix4d transform = Matrix4d::Identity();
Vector3d point{1.0, 2.0, 3.0};
Vector3d result = transform.topLeftCorner<3, 3>() * point;

// Dynamic-size — heap-allocated, still vectorized
MatrixXd A = MatrixXd::Random(100, 100);
MatrixXd B = MatrixXd::Random(100, 50);
MatrixXd C = A * B;  // single expression, no temporary matrices

// Solvers
Eigen::ColPivHouseholderQR<MatrixXd> qr(A);
VectorXd x = qr.solve(b);  // Ax = b

// Eigenvalues
Eigen::SelfAdjointEigenSolver<MatrixXd> solver(A);
auto eigenvalues = solver.eigenvalues();

Expression templates — zero-cost abstractions

cpp

Godbolt

// This expression:
MatrixXd result = A * B + C;

// Does NOT create A*B as a temporary, then add C.
// Eigen generates a single fused loop:
// result(i,j) = sum_k(A(i,k)*B(k,j)) + C(i,j)
// This is critical for performance — eliminates O(n^2) allocations

BLAS / LAPACK

For maximum performance on CPU, use optimized BLAS implementations:

Implementation	Vendor	Best for
OpenBLAS	Community	AMD/Intel, open source
Intel MKL / oneMKL	Intel	Intel CPUs (especially Xeon)
BLIS	UT Austin	Research, modular
Apple Accelerate	Apple	macOS/iOS (free, well-tuned)

cmake

# CMake — find BLAS
find_package(BLAS REQUIRED)
find_package(LAPACK REQUIRED)
target_link_libraries(myapp PRIVATE BLAS::BLAS LAPACK::LAPACK)

# Or: link Intel MKL explicitly
target_link_libraries(myapp PRIVATE
    -Wl,--start-group
    mkl_intel_lp64 mkl_sequential mkl_core
    -Wl,--end-group -lm -lpthread)

cpp

Godbolt

// Direct BLAS calls (when Eigen abstraction insufficient)
#include <cblas.h>

// DGEMM: C = alpha * A*B + beta * C
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
    m, n, k,
    alpha, A, lda,
    B, ldb,
    beta, C, ldc);

OpenMP — shared-memory parallelism

cpp

Godbolt

#include <omp.h>
#include <vector>
#include <cmath>

// Parallel for loop — divide work across CPU cores
void compute_sqrt(std::vector<double>& data) {
    #pragma omp parallel for schedule(static)
    for (std::size_t i = 0; i < data.size(); ++i) {
        data[i] = std::sqrt(data[i]);
    }
}

// Reduction — safe summation
double dot_product(const double* a, const double* b, int n) {
    double sum = 0.0;
    #pragma omp parallel for reduction(+:sum)
    for (int i = 0; i < n; ++i)
        sum += a[i] * b[i];
    return sum;
}

// Task parallelism
#pragma omp parallel
{
    #pragma omp single
    {
        #pragma omp task
        compute_subtree(left);
        #pragma omp task
        compute_subtree(right);
        #pragma omp taskwait
    }
}

cmake

find_package(OpenMP REQUIRED)
target_link_libraries(myapp PRIVATE OpenMP::OpenMP_CXX)

MPI — distributed-memory parallelism

MPI enables computation across multiple nodes in a cluster.

cpp

Godbolt

#include <mpi.h>
#include <vector>

int main(int argc, char** argv) {
    MPI_Init(&argc, &argv);

    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);   // this process's ID
    MPI_Comm_size(MPI_COMM_WORLD, &size);   // total number of processes

    // Distribute work — each rank processes a slice
    int n = 1000000;
    int local_n = n / size;
    std::vector<double> local_data(local_n);

    // Fill local data...
    for (int i = 0; i < local_n; ++i)
        local_data[i] = rank * local_n + i;

    // Reduce: sum all local_sums into global_sum on rank 0
    double local_sum = 0;
    for (double x : local_data) local_sum += x;

    double global_sum = 0;
    MPI_Reduce(&local_sum, &global_sum, 1, MPI_DOUBLE,
               MPI_SUM, 0, MPI_COMM_WORLD);

    if (rank == 0)
        std::println("Total sum: {}", global_sum);

    MPI_Finalize();
    return 0;
}

bash

mpic++ -O3 -std=c++20 -o mpi_app main.cpp
mpirun -np 8 ./mpi_app   # run on 8 processes

CUDA C++ — GPU computing

cpp

Godbolt

#include <cuda_runtime.h>

// Kernel — runs on GPU
__global__ void vector_add(const float* a, const float* b, float* c, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n)
        c[idx] = a[idx] + b[idx];
}

int main() {
    int n = 1 << 20;  // 1M elements
    std::size_t bytes = n * sizeof(float);

    // Allocate unified memory (accessible on CPU and GPU)
    float *a, *b, *c;
    cudaMallocManaged(&a, bytes);
    cudaMallocManaged(&b, bytes);
    cudaMallocManaged(&c, bytes);

    // Init on CPU
    for (int i = 0; i < n; ++i) { a[i] = i; b[i] = i * 2; }

    // Launch kernel — 256 threads per block
    int threads = 256;
    int blocks = (n + threads - 1) / threads;
    vector_add<<<blocks, threads>>>(a, b, c, n);

    cudaDeviceSynchronize();  // wait for GPU

    cudaFree(a); cudaFree(b); cudaFree(c);
}

bash

nvcc -O3 -std=c++17 -o cuda_app main.cu

SIMD intrinsics

For maximum throughput on specific operations:

cpp

Godbolt

#include <immintrin.h>  // AVX2

// Multiply two arrays of 8 floats at once (AVX)
void multiply_avx(const float* a, const float* b, float* c, int n) {
    int i = 0;
    for (; i + 8 <= n; i += 8) {
        __m256 va = _mm256_loadu_ps(a + i);
        __m256 vb = _mm256_loadu_ps(b + i);
        _mm256_storeu_ps(c + i, _mm256_mul_ps(va, vb));
    }
    // Scalar tail
    for (; i < n; ++i) c[i] = a[i] * b[i];
}

Better: use std::experimental::simd (C++23) or Eigen's auto-vectorization and avoid manual intrinsics.

Memory layout — row-major vs column-major

This is the #1 performance mistake in scientific code:

cpp

Godbolt

// Row-major (C++ default) — iterate rows in inner loop
// Cache-friendly for A[i][j] where j varies fastest
double A[N][N];
for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j)
        A[i][j] = compute(i, j);  // sequential — cache friendly

// BLAS uses column-major (Fortran convention)
// Eigen defaults to column-major internally but handles correctly
// std::mdspan (C++23) explicitly models both
#include <mdspan>

auto md = std::mdspan(data, std::extents<int, N, N>{});
// Column-major:
auto cm = std::mdspan(data,
    std::layout_left::mapping<std::extents<int, N, N>>{});

Scientific library ecosystem

Library	Purpose
Eigen	Dense linear algebra, geometry
Armadillo	MATLAB-like syntax over BLAS
blaze	High-performance dense/sparse math
GSL	GNU Scientific Library — vast algorithm collection
FFTW	Fastest Fourier Transform in the West
Boost.Multiprecision	Arbitrary precision arithmetic
deal.II	Finite element analysis
VTK	Scientific visualization
HDF5	Parallel I/O for large datasets
NetCDF	Array-oriented scientific data
Kokkos	Performance-portable parallel programming
RAJA	Portability layer (OpenMP/CUDA/HIP)

Kokkos — performance portability

Kokkos lets you write parallel code once that runs correctly on CPU (OpenMP threads), NVIDIA GPU (CUDA), AMD GPU (HIP), and Intel GPU (SYCL) by selecting a backend at compile time. The key abstractions are execution spaces (where code runs) and memory spaces (where data lives).

cpp

Godbolt

#include <Kokkos_Core.hpp>

int main(int argc, char** argv) {
    Kokkos::initialize(argc, argv);
    {
        const int N = 1'000'000;

        // View<double*> — a managed N-element array in the default memory space
        // (host RAM for OpenMP backend, device VRAM for CUDA backend)
        Kokkos::View<double*> a("a", N);
        Kokkos::View<double*> b("b", N);
        Kokkos::View<double*> c("c", N);

        // Initialize on device — same code runs as CUDA kernels or OpenMP threads
        Kokkos::parallel_for("init", N, KOKKOS_LAMBDA(int i) {
            a(i) = i * 1.0;
            b(i) = i * 2.0;
        });

        // Element-wise addition
        Kokkos::parallel_for("add", N, KOKKOS_LAMBDA(int i) {
            c(i) = a(i) + b(i);
        });

        // Reduction — sum all elements of c
        double sum = 0.0;
        Kokkos::parallel_reduce("sum", N,
            KOKKOS_LAMBDA(int i, double& local_sum) {
                local_sum += c(i);
            }, sum);
        // parallel_for/reduce are asynchronous — fence before using result on CPU
        Kokkos::fence();

        std::println("sum = {}", sum);
    }
    Kokkos::finalize();
}

2D View and matrix operations

View supports multi-dimensional layouts. LayoutLeft is column-major (Fortran/BLAS order), LayoutRight is row-major (C order, cache-friendly for row iteration).

cpp

Godbolt

const int M = 1000, N = 1000;

// Column-major 2D view (good for BLAS interop)
Kokkos::View<double**, Kokkos::LayoutLeft> A("A", M, N);
Kokkos::View<double*>                      x("x", N);
Kokkos::View<double*>                      y("y", M);

// Matrix-vector multiply: y = A * x
Kokkos::parallel_for("matvec", M, KOKKOS_LAMBDA(int i) {
    double sum = 0.0;
    for (int j = 0; j < N; ++j)
        sum += A(i, j) * x(j);
    y(i) = sum;
});

Explicit host/device copies

When a View must be on the device but you need to read results on the host, use a mirror:

cpp

Godbolt

// Allocate device View
Kokkos::View<double*> dev("device_data", N);

// Create a host mirror (no-op if already on host)
auto host = Kokkos::create_mirror_view(dev);

// Fill on host, copy to device
for (int i = 0; i < N; ++i) host(i) = i * 3.14;
Kokkos::deep_copy(dev, host);  // host -> device transfer

// ... run kernel on dev ...

// Copy results back to host for inspection
Kokkos::deep_copy(host, dev);
Kokkos::fence();
std::println("result[0] = {}", host(0));

CMake integration

cmake

# Clone or use a system install of Kokkos
find_package(Kokkos REQUIRED)
target_link_libraries(myapp PRIVATE Kokkos::kokkos)

# Configure backend at cmake time:
# cmake -DKokkos_ENABLE_CUDA=ON -DKokkos_ARCH_AMPERE80=ON ..  # NVIDIA A100
# cmake -DKokkos_ENABLE_OPENMP=ON ..                           # CPU threads
# cmake -DKokkos_ENABLE_HIP=ON ..                              # AMD GPU

C++17 parallel algorithms

The standard library's parallel algorithms (<algorithm>, <numeric>) accept an execution policy. With a supported library (Intel TBB, NVHPC stdpar), std::execution::par_unseq vectorizes and parallelizes on CPU cores or GPUs — no OpenMP pragmas needed.

cpp

Godbolt

#include <algorithm>
#include <execution>
#include <numeric>
#include <vector>

std::vector<double> data(10'000'000);
std::iota(data.begin(), data.end(), 0.0);

// Parallel sort — automatically uses thread pool
std::sort(std::execution::par_unseq, data.begin(), data.end());

// Parallel transform: apply sqrt to every element in-place
std::transform(std::execution::par_unseq,
               data.begin(), data.end(), data.begin(),
               [](double x) { return std::sqrt(x); });

// Parallel reduction — sum
double total = std::reduce(std::execution::par_unseq,
                           data.begin(), data.end(), 0.0);

// Parallel transform_reduce — dot product
std::vector<double> a(N), b(N);
double dot = std::transform_reduce(std::execution::par_unseq,
    a.begin(), a.end(), b.begin(), 0.0,
    std::plus<>{}, std::multiplies<>{});

cmake

# Link TBB for CPU parallelism on GCC/Clang:
find_package(TBB REQUIRED)
target_link_libraries(myapp PRIVATE TBB::tbb)

Execution policies:

std::execution::seq — sequential (default behavior)
std::execution::par — parallel, no vectorization
std::execution::par_unseq — parallel + SIMD vectorization (requires the lambda to be safe to vectorize: no locks, no memory aliasing)

FFTW — fast Fourier transforms

FFTW (Fastest Fourier Transform in the West) achieves near-peak FLOP throughput by benchmarking multiple algorithms at plan creation time and choosing the fastest for your specific hardware and array size.

cpp

Godbolt

#include <fftw3.h>
#include <complex>
#include <vector>

int N = 1024;

// Allocate FFTW-aligned memory (required for SIMD)
auto* in  = fftw_alloc_complex(N);
auto* out = fftw_alloc_complex(N);

// Plan: FFTW_MEASURE benchmarks internally to pick the fastest algorithm
fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_MEASURE);

// Fill input
for (int i = 0; i < N; ++i) {
    in[i][0] = std::cos(2 * M_PI * i / N);  // real part
    in[i][1] = 0.0;                          // imaginary part
}

// Execute — can be called repeatedly with the same plan
fftw_execute(plan);

// out[k][0] + i*out[k][1] is the k-th frequency component
double magnitude_dc = std::hypot(out[0][0], out[0][1]);

fftw_destroy_plan(plan);
fftw_free(in);
fftw_free(out);

cmake

find_package(PkgConfig REQUIRED)
pkg_check_modules(FFTW REQUIRED fftw3)
target_link_libraries(myapp PRIVATE ${FFTW_LIBRARIES})
target_include_directories(myapp PRIVATE ${FFTW_INCLUDE_DIRS})

Tips:

Create plans once; reuse for repeated transforms of the same size
FFTW_ESTIMATE skips benchmarking (fast plan creation, slower transforms) — useful in tests
FFTW_MEASURE or FFTW_PATIENT benchmark at startup — best for production
FFTW supports multi-threaded transforms via fftw_init_threads() + fftw_plan_with_nthreads(n)

Profiling HPC code

Identifying bottlenecks before optimizing is essential — intuition is unreliable for memory-bound or vectorization-limited code.

bash

# Linux perf — CPU cycles, cache misses, branch mispredictions
perf stat ./myapp
perf record -g ./myapp && perf report   # flame graph

# Valgrind cachegrind — cache miss analysis
valgrind --tool=cachegrind ./myapp
cg_annotate cachegrind.out.*            # per-line cache stats

# NVIDIA GPU profiling
nsys profile --stats=true ./myapp       # Nsight Systems — timeline
ncu --set full ./myapp                   # Nsight Compute — kernel roofline

# Intel VTune (CPU + GPU profiling on Intel hardware)
vtune -collect hotspots ./myapp
vtune -collect memory-access ./myapp    # memory bandwidth analysis

# OpenMP thread timeline
OMP_NUM_THREADS=8 KMP_AFFINITY=granularity=fine,compact ./myapp

Roofline model — is your kernel compute- or memory-bound?

The roofline model plots arithmetic intensity (FLOP/byte) against peak FLOP/s and peak memory bandwidth to show which limit applies:

cpp

Godbolt

Peak FLOP/s (e.g. 10 TFLOP/s)  ──────────────────────────────────── compute bound
                                                              /
                                                             /
                                                            /  memory bound
Peak BW × intensity  ──────────────────────────────────── /
                      0.1          1            10         AI (FLOP/byte)

If your kernel sits left of the ridge point (memory-bound), optimize data layout and access patterns. If right (compute-bound), look at instruction throughput and vectorization.