We use Amdahl's Law (Summary Section 1.3.2): $SU(p) = \frac{1}{(1-f) + \frac{f}{p}}$.

$4 = \frac{1}{(1-f) + \frac{f}{16}}$

$(1-f) + \frac{f}{16} = \frac{1}{4} = 0.25$

$1 - 0.25 = f - \frac{f}{16}$

$0.75 = f(1 - \frac{1}{16}) = \frac{15f}{16}$

$f = \frac{0.75 \times 16}{15} = \frac{12}{15} = 0.8$

80% of the execution time is parallelized.

The maximum speedup (Amdahl's Limit) occurs as $p \to \infty$.

$SU(\infty) = \frac{1}{1-f}$ (Section 1.3.2).

$SU(\infty) = \frac{1}{1-0.8} = \frac{1}{0.2} = 5$.

The expected speedup is 5.

(Section 6.3.1)

1. Synchronization Overhead: Time spent coordinating threads, such as waiting at barriers or acquiring locks (Section 6.3.2).
2. Load Imbalance: Uneven distribution of work, causing some processing units to remain idle while waiting for others to finish.

(Section 1.5)

• Data Parallelism (1.5.1): The same computation is performed concurrently on different subsets of the data (Split data, same computation).
• Functional Parallelism (1.5.2): Different computations (tasks) are performed concurrently (Split computation).

A typical HPC architecture involves the following hierarchy:

• Physical Structure: A Server Rack houses multiple Nodes (the blades). Nodes require infrastructure like a Power supply and Water cooling.
• Node Internals (Mainboard): The mainboard features one or more NUMA Sockets, each holding a Processor (CPU). Main memory (DRAM) is attached to the sockets.
• Processor Internals: The processor contains cores and Cache.
• Execution: Cores execute Hardware Threads. The operating system runs a Process, which consists of Software Threads. In MPI, a process is identified by its Rank.

In modern multi-socket systems, memory is physically distributed, with memory banks attached locally to specific sockets. Accessing local memory is fast (low latency). Accessing memory attached to a remote socket requires traversing the interconnect between the sockets, which incurs significantly higher latency.

This effect is called NUMA (Non-Uniform Memory Access) (Section 1.6.2.1).

False sharing occurs when different cores modify independent variables located on the same cache line, causing excessive coherence traffic between their private caches (Section 6.4). Core 0 and Core 2 have private L1 (16B line) and L2 (32B line) caches. To prevent false sharing, we must align the data based on the largest private cache line size involved, which is 32 bytes.

uint32_t is 4 bytes. Padding needed = 32 bytes (Line Size) - 4 bytes (variable size) = 28 bytes.

struct CounterData {
    uint32_t counterThreadZero;

    // Padding to ensure variables are on separate 32-byte cache lines
    char padding[28];

    uint32_t counterThreadTwo;
};

System context: 2 NUMA sockets, 4 cores each (8 total). 8 threads used.

1. NUMA Effects / First Touch Policy (Section 6.5.2): The `image` data is initialized sequentially by the main thread (Lines 35-39). Due to the First Touch policy, all memory pages are physically allocated on the NUMA node where the main thread runs. During the parallel computation (Lines 44-46), threads running on the remote socket will experience high latency due to remote memory access.
   • Solution: Initialize the `image` data in parallel using the same partitioning scheme used for the computation.
2. Poor Spatial Locality / Strided Access (Section 6.4.1): In `normalize_image` (Lines 16-21), the loops are nested such that the inner loop iterates over `row`. Since the image is stored in row-major order (`image[row * IMG_WIDTH + col]`), the inner loop accesses memory with a large stride (`IMG_WIDTH`). This leads to poor cache utilization.
   • Solution: Interchange the loops so the inner loop iterates over `col`, ensuring contiguous memory access.

• Advantage: Low overhead for creation and context switching, as they do not require expensive system calls (kernel involvement).
• Disadvantage: A blocking system call made by one user-level thread may block the entire process, as the OS scheduler is typically unaware of individual user-level threads.
• Use Case: Implementing fine-grained concurrency models like fibers or coroutines, often used in event-driven applications.

(Section 2.4.2)

• Spin-Lock: A thread waits by continuously polling the lock variable ("busy-waiting").
  • Advantage: Low latency (fast acquisition when the lock is released).
  • Disadvantage: Wastes CPU cycles while spinning.
• Yielding Lock: A thread yields its execution slice to the OS scheduler (goes to sleep) if the lock is occupied.
  • Advantage: Does not waste CPU resources while waiting.
  • Disadvantage: High latency due to the overhead of context switching and rescheduling.

The problem is a Loop-Carried Flow Dependence (Section 4.2.1). Statement `c[i] = a[i-1] + b[i-1]` at iteration $i$ reads `a[i-1]`, which is written by the statement `a[i] = a[i] * b[i]` in the previous iteration $i-1$. This dependency across iterations prevents the loop from being vectorized, as all iterations in a vector operation must be independent.

We use loop alignment (Section 4.2.2) to shift the execution so the dependency occurs within the same iteration (loop-independent). We align the computation of `a[i]` with `c[i+1]`. This requires peeling the first iteration for the `c` calculation and the last iteration for the `a` calculation.

void operation_aligned(float *a, float *b, float *c, int size) {
    if (size <= 1) return;

    // Peel first iteration for the 'c' calculation (original i=1)
    c[1] = a[0] + b[0];

    // Aligned loop (runs from i=1 to size-2)
    for(int i = 1; i < size - 1; i++){
        a[i] = a[i] * b[i];       // S1(i)
        // Now uses the value of a[i] computed in this iteration
        c[i+1] = a[i] + b[i];     // S2(i+1)
    }

    // Peel last iteration for the 'a' calculation (original i=size-1)
    a[size-1] = a[size-1] * b[size-1];
}

We vectorize the aligned loop. VL = 8 (256-bit / 32-bit float). We use unaligned intrinsics.

#include 

void operation_vectorized(float *a, float *b, float *c, int size){
    if (size <= 1) return;

    // Peeling (scalar)
    c[1] = a[0] + b[0];

    int i = 1;
    const int VL = 8;

    // Main Vector Loop (from i=1 up to the last full vector in the size-2 range)
    for(; i <= (size - 2) - VL; i += VL){
        __m256 a256, b256, c256;

        // Load a[i..i+7] and b[i..i+7] (Unaligned)
        a256 = _mm256_loadu_ps(&a[i]);
        b256 = _mm256_loadu_ps(&b[i]);

        // S1(i): a = a * b
        a256 = _mm256_mul_ps(a256, b256);

        // Store updated a
        _mm256_storeu_ps(&a[i], a256);

        // S2(i+1): c = a + b (using the newly computed a)
        c256 = _mm256_add_ps(a256, b256);

        // Store c[i+1..i+8]
        _mm256_storeu_ps(&c[i+1], c256);
    }

    // Remainder Loop (Scalar)
    for( ; i < size - 1; i++ ){
        a[i] = a[i] * b[i];
        c[i+1] = a[i] + b[i];
    }

    // Peeling (scalar)
    a[size-1] = a[size-1] * b[size-1];
}

The workload per row in a sparse matrix is irregular, leading to potential load imbalance. We must use `dynamic` scheduling (Section 3.3.3) for efficiency. We also need a reduction for the `max` variable (Section 3.6) and ensure loop-local variables (`nnz`, `num_nnz`) are private.

#pragma omp parallel for reduction(max:max) private(nnz, num_nnz) schedule(dynamic)

We implement a manual reduction using thread-private variables and a critical section for the final aggregation (Section 3.5.1).

int calculate_sum2(double *data, int size){
    double sum = 0.0;
    /* OMP PRAGMA1 */
    #pragma omp parallel
    {
        // Note: Template used 'int local_sum', corrected to 'double' as data is double.
        double local_sum = 0.0;
        /* OMP PRAGMA2 */
        #pragma omp for
        for (int i = 0; i < size; i++) {
            local_sum += data[i];
        }
        /* OMP PRAGMA3 */
        #pragma omp critical
        {
            sum += local_sum;
        }
    }
    return sum;
}

(Section 3.1.2)

1. Compiler Directives (Pragmas, e.g., #pragma omp parallel).
2. Runtime Library Routines (Functions, e.g., omp_get_thread_num()).
3. Environment Variables (e.g., OMP_NUM_THREADS).

Size=12 (Ranks 0-11). `col = rank / 4`; `row = rank % 4`.

We analyze Rank 6: `col`=6/4=1. `row`=6%4=2.

Communicator c1: `MPI_Comm_split(color=row, key=col)` (Section 9.8.2). Rank 6 is in group `color=2` (all ranks where `rank % 4 == 2`). Members: {2, 6, 10}. The keys for these members are their `col` values: Rank 2 (col=0), Rank 6 (col=1), Rank 10 (col=2). Ranks in the new communicator are assigned based on ascending key order: Rank 2 gets new rank 0, Rank 6 gets new rank 1, Rank 10 gets new rank 2. So, `rank_c1` for Rank 6 is 1. `Allgather` collects the world ranks {2, 6, 10}.

Communicator c2: `MPI_Comm_split(color=col, key=4-row)`. Rank 6 is in group `color=1` (all ranks where `rank / 4 == 1`). Members: {4, 5, 6, 7}. The keys for these members are `4-row`: Rank 4 (4-0=4), Rank 5 (4-1=3), Rank 6 (4-2=2), Rank 7 (4-3=1). Ranks are assigned based on ascending key order: Rank 7 (key=1) gets new rank 0, Rank 6 (key=2) gets new rank 1, Rank 5 (key=3) gets new rank 2, Rank 4 (key=4) gets new rank 3. So, `rank_c2` for Rank 6 is 1. `Allgather` collects the world ranks in the new order: {7, 6, 5, 4}.

Output (at Rank 6):

rank_c1: 1 members: 2 6 10 
rank_c2: 1 members: 7 6 5 4 

• Blocking routines (e.g., MPI_Send): Return only after the operation is locally complete, meaning the communication buffer is safe to reuse (Section 8.4.2).
• Nonblocking routines (e.g., MPI_Isend): Return immediately after initiating the operation. The application must later verify completion (e.g., using MPI_Wait) before reusing the buffer (Section 9.5).

#include 
// assume print() is defined

int main(int argc, char *argv[])
{
    MPI_Init(&argc, &argv);
    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    int local_results[] = { rank, rank * 2, rank * 2}; // 3 elements
    
    // TODO: declare recvbuf
    // Needs to hold the result (3 integers), only significant on root
    int recvbuf[3];

    // TODO: reduce local_results into recvbuf using the sum operation at rank 0
    // MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm)
    MPI_Reduce(local_results, recvbuf, 3, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);

    if (rank == 0) {
        // print content of recvbuf
        // print(recvbuf);
    }
    MPI_Finalize();
    return 0;
}

Ranks 0, 1, 2, 3. `local_results = {rank, rank*2, rank*2}`.

• Rank 0: {0, 0, 0}
• Rank 1: {1, 2, 2}
• Rank 2: {2, 4, 4}
• Rank 3: {3, 6, 6}

The sum is performed element-wise:

Sum: {0+1+2+3, 0+2+4+6, 0+2+4+6} = {6, 12, 12}.

The content of `recvbuf` at rank 0 is {6, 12, 12}.

We must use `MPI_Init_thread` to initialize MPI with thread support (Section 9.13.1). For the typical hybrid model where MPI calls are made outside OpenMP regions by the main thread, `MPI_THREAD_FUNNELED` is sufficient and appropriate (Section 9.13.4).

#include 
#include 

// Forward declaration of the function and Pixel struct
struct Pixel; 
void average_luminance(Pixel *pixels, int size, int root);

int main(int argc, char *argv[]) {
    int required = MPI_THREAD_FUNNELED;
    int provided;

    MPI_Init_thread(&argc, &argv, required, &provided);

    if (provided < required) {
        fprintf(stderr, "Error: MPI implementation does not provide required thread support.\n");
        MPI_Abort(MPI_COMM_WORLD, 1);
    }

    // Assume pixels, size, root are initialized here for the example
    // average_luminance(pixels, size, root);

    MPI_Finalize();
    return 0;
}

We use a hybrid approach: MPI Scatter to distribute the data, OpenMP parallel for reduction to calculate local sums, and MPI Reduce to aggregate the global sum at the root.

#include 
#include 
#include 
#include 

struct Pixel {
    unsigned short r;
    unsigned short g;
    unsigned short b;
};

void average_luminance(Pixel *pixels, int size, int root){
    int commSize, commRank;
    MPI_Comm_size(MPI_COMM_WORLD, &commSize);
    MPI_Comm_rank(MPI_COMM_WORLD, &commRank);

    // 1. Define MPI Datatype for Pixel (3 contiguous unsigned shorts)
    MPI_Datatype MPI_PIXEL;
    MPI_Type_contiguous(3, MPI_UNSIGNED_SHORT, &MPI_PIXEL);
    MPI_Type_commit(&MPI_PIXEL);

    int local_size = size / commSize;

    // 2. Allocate local buffer
    Pixel *localPixels = (Pixel*) malloc(local_size * sizeof(Pixel));

    // 3. Distribute data from root to all processes
    MPI_Scatter(pixels, local_size, MPI_PIXEL, localPixels, local_size, MPI_PIXEL, root, MPI_COMM_WORLD);

    // 4. Local computation (OpenMP Parallelization)
    double local_sum = 0.0;
    #pragma omp parallel for reduction(+:local_sum) schedule(static)
    for(int i = 0; i < local_size; i++){
        double luminance = 0.2126 * localPixels[i].r
                         + 0.7152 * localPixels[i].g
                         + 0.0722 * localPixels[i].b;
        local_sum += luminance;
    }

    // 5. Global reduction to sum up local_sum from all ranks at root
    double global_sum;
    MPI_Reduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, root, MPI_COMM_WORLD);

    // 6. Output at root
    if (commRank == root) {
        printf("Average luminance: %f\n", global_sum / size);
    }

    // Cleanup
    free(localPixels);
    MPI_Type_free(&MPI_PIXEL);
}

Hardware Context: 2 sockets, 4 cores/socket, 2-way SMT = 16 hardware threads/node.
Configuration: 4 MPI processes per node.

Performance Issues:

1. Massive OpenMP Oversubscription: `OMP_NUM_THREADS = '16'`. Each of the 4 MPI processes attempts to launch 16 threads, resulting in 64 total software threads competing for 16 hardware threads. This causes excessive context switching overhead.
2. Lack of Thread Pinning (Poor Locality): `OMP_PROC_BIND = 'false'`. Threads are not bound to cores (Section 6.5.3). They can migrate freely across the node, including across NUMA boundaries, which destroys data locality and cache performance.
3. Suboptimal Hybrid Strategy: Running 4 MPI processes on a 2-socket node (2 per socket) is generally less efficient than the recommended hybrid strategy of 1 MPI process per socket (Section 9.13.4), as it leads to resource contention for memory bandwidth within each socket.

Fixes (Optimal Strategy):

The best approach is to reconfigure the execution to match the hardware: 1 MPI process per socket (2 total per node), and have OpenMP utilize the cores within each socket. We prioritize physical cores (4 per socket) as SMT is often not beneficial for HPC.

1. Fix MPI Execution Command (assuming OpenMPI):

# Start 2 processes per node, mapping one to each socket
mpirun --map-by socket --bind-to socket -np 2 ./your_program

2. Set OpenMP Environment Variables:

export OMP_NUM_THREADS=4       # Utilize the 4 physical cores per socket
export OMP_PROC_BIND=close     # Keep threads close within the socket's cores
export OMP_PLACES=cores        # Define the places for binding as cores