Parallel Computing with MPI — The Sieve of Eratosthenes and Performance Optimization

Overview

Target:

• Implementing the parallel Eratosthenes Sieve algorithm using MPI programming.

• Analyze and optimize program performance. There are three basic optimizations:

  • optimize1: eliminate even numbers
  • optimize2: eliminate broadcast
  • optimize3: improve cache utilization

• Finally, optimize4 is for other optimization tricks. Including using memset, unroll loops, __builtin_popcountll, register keyword, inline assembly, etc.

Experimental Environment:

Huawei Cloud Elastic Server, specifications/image:

_Kunpeng General Computing-Plus | 2vCPUs | 4GiB | kc1.large.2_ > _Ubuntu 18.04 server 64bit with ARM(40GiB)_

Optimizations

$\boxed{\text{base}}$：Normal MPI Eratosthenes Sieve

Before discussing parallelization and various optimizations, we need first understand the implementation principle of the Eratosthenes Sieve:

• For any positive integer $n\gt 1$, its multiple $nx$ is composite when $x\gt 1$. Consider each number $a$ in ascending order and mark its multiples $ax$ as composites. The unmarked numbers are prime.

The serial sieve is straightforward. Its MPI parallelization follows naturally.

• With $p$ processes, divide the $n$ numbers into $p$ segments. Each process handles only its assigned interval. More specifically, process $i$ handles the interval from $\lfloor in/p\rfloor$ to $\lfloor(i+1)n/p\rfloor-1$.

  Since 0 and 1 do not need to be processed, the actual data starts at 2. In the code, the bounds are low_value = 2 + i * (n - 1) / p; and high_value = 1 + (id + 1) * (n - 1) / p;.

• What each process need to do is taking a previously identified prime $prime$, finds the first multiple of $prime$ in its interval and records its index as $first$. Then it repeatedly adds $prime$ to $first$ to mark composites. After finishing the interval, it moves to the next prime.

• The sieving primes come from broadcasts by process 0. Process 0 identifies each new prime and broadcasts it to all other processes.

• The marking array only need to be as large as the local interval. Map $first$ to this array as follows:

  • If $prime^2\gt low_{value}$, where $low_{value}$ is the lower bound of the process's interval, use $first = prime^2 - low_{value}$

    Why start at $prime^2$? Let the preceding multiple be $prime\cdot a$. Then we have $a\lt prime$.

    • If $a$ is prime, it has already been used for sieving because the primes are processed in ascending order. Therefore $prime\cdot a$ has already been marked.
    • If $a$ is not prime, all its prime factors are smaller than $prime$. At least one of those previously processed primes has already marked $prime\cdot a$

  • If $prime^2 \le low_{value}$:

    • If $low_{value}\mod prime \neq 0$ $first = prime-(low_{value}\mod prime)$
    • Otherwise, let $first=0$ 。

Flowchart for base.cpp:

Result for base.cpp:

As the number of processes increases, the reduction in runtime becomes smaller. Communication efficiency may be a factor. Amdahl's law also limits the speedup according to how much of the work can be parallelized. The parallel portion of base.cpp still has considerable room for optimization.

$\boxed{\text{optimize1}}$: Eliminate Even Numbers

Obviously, even number other than $2$ are not prime, so we can eliminate them in advance and halve the problem size. Compared with base.cpp:

• Each interval is half as large, and the marked array requires half the storage.

• Calculating $first$ requires some additional care because the array has been halved. The overall approach is unchanged:

  if (prime * prime > low_value)
      first = (prime * prime - low_value) / 2; // Halve the index
  else {
      if (!(low_value % prime))
          first = 0;
      else if (low_value % prime % 2 == 0)
          first = prime - ((low_value % prime) / 2);
      else
          first = (prime - (low_value % prime)) / 2;
  }
  

  The case where $low_{value}\ \%\ prime$ is even deserves an explanation:

  • When the remainder is even, $prime\cdot a-low_{value}$ is odd, where $prime\cdot a$ is the first multiple of prime greater than low_value. Since even candidates have been excluded, low_value must be odd, which makes $prime\cdot a$ even. We must skip that multiple. Advancing to the next odd multiple gives an odd composite that needs to be marked; successive odd multiples are separated by $2\cdot prime$ in value.

    

Flowchart for optimize1.cpp

The structure is the same as the base.cpp flowchart. The initial value of prime is 3 rather than 2.

Result for optimize1.cpp

The runtime is approximately half that of base.cpp, showing that removing even numbers is a direct and effective improvement.

$\boxed{\text{optimize2}}$: Eliminate Broadcasts

Each process currently obtains its sieving primes from process $0$ through broadcasts, but these broadcasts can be eliminated entirely.

Each process can generate the primes it needs independently. Sieving up to $n$ only requires primes up to $\sqrt n$. Generating those primes with a serial algorithm removes the need to wait for broadcasts.

I used Euler's sieve to generate these primes because it is simple to implement and has linear time complexity, $O(n)$ for a sieve limit of $n$. Its key rule is to stop the inner loop when $i\ \%\ prime=0$ while marking $i\cdot prime$.

• The primes are processed in ascending order. Once a prime divides $i$, products involving subsequent primes will be marked elsewhere using that smaller prime. For example, with $i=2\times3=6$ and $prime=2$, the next prime is $3$. The product $i\cdot3=18$ will already be marked when $i=9$ and $prime=2$.

      vector<int> pri;
      int pri_i = 1;
      int sqrt_N = sqrt(n) + 10;
      bool not_prime[sqrt_N + 10] = {0};
      for (int i = 2; i <= sqrt_N; i++) {
          if (!not_prime[i]) pri.push_back(i);
          for (int prime_i : pri) {
              if (i * prime_i > sqrt_N) break;
              not_prime[i * prime_i] = true;
              if (i % prime_i == 0) break; //Break immediately
          }
      }
  

Flowchart for optimize2.cpp:

Generate the primes up to $\sqrt n$ beforehand. Each process then obtains the next prime from its own list instead of process $0$.

Results for optimize2.cpp:

The main difference from optimize1.cpp appears at 16 processes. This suggests that removing communication becomes more beneficial as the number of processes grows.

$\boxed{\text{optimize3}}$: Cache optimization

Cache is small, fast memory. Accessing cached data is much faster than accessing main memory. Processing data within a working set that fits in cache can improve the hit rate and overall performance. For an array, accessing an element brings a cache line containing nearby elements into the cache for subsequent use.

Nearby sieving primes also often mark the same composite within a given interval. Combining this observation with cache locality suggests splitting each process's data into smaller blocks based on the cache size.

For op2, each process continuously uses primes $2,3,5...$ to sieve the whole interval. For op3, we divide the interval into blocks. Each process first sieve a block and then continue to the next block.

I chose the block size as follows:

• First, inspect the cache size using lscpu. The L3 cache size is $32768\times1024=33554432\ B$.

• The final evaluation uses 8 processes, so first divide the cache evenly: $33554432/8=4194304\ B$ per process. To leave headroom, use one quarter of that amount, or $1048576\ B$.

  The implementation uses a char array for marks and a long long int array for sieving primes, with primes needed up to $\sqrt n$. As a rough estimate, $x$ char entries are associated with $\sqrt x$ long long int entries. I approximated the average storage per candidate as $1\ (\text{sizeof char})+\sqrt{8\ (\text{sizeof long long int})}\approx4\ B$.

  This gives $Block_{size}=1048576/4=262144$.

  In code, the resulting rule is CACHE_SIZE / p / 4 / 4. I tried several formulas and sizes; this one appeared to work well, so I kept it.

• Processing each block is similar to the earlier approach: subdivide the interval again, calculate the block bounds, sieve, and count the remaining primes.

  • Number of candidates: N = (n - 1) / 2.

    Block quantity as stated in the source (indexed from 0): Block_N = CacheBlock_size - 1.

    Index of the first candidate assigned to this process: low_index = id * (N / p) + MIN(id, N % p).

    Index of the last candidate assigned to this process: high_index = (id + 1) * (N / p) + MIN(id + 1, N % p) - 1.

    Process lower bound: low_value = low_index * 2 + 3.

    Process upper bound: high_value = high_index * 2 + 3.

    Index of the first candidate in the block: Block_first = Block_id * Block_N + MIN(Block_id, Block_rest) + low_index.

    Index of the last candidate in the block: Block_last = (Block_id + 1) * Block_N + MIN(Block_id + 1, Block_rest) - 1 + low_index.

    Block lower bound: Block_low_value = 2 * Block_first + 3.

    Block upper bound: Block_high_value = 2 * Block_last + 3, except for the final block, where Block_high_value = high_value.

  • Continue as before, leaving the inner loop when $prime^2>Block_{high\ value}$.

Flowchart for optimize3.cpp:

Results for optimize3.cpp:

Compared with optimize2.cpp, the speedup is substantial. Efficient cache access makes a clear difference.

$\boxed{\text{optimize4}}$: Final Optimizations

This stage allowed more experimentation. I describe the optimization process roughly in chronological order. Most development tests below used 16 processes and a problem size of $3\times10^9$.

1. Initialize marked with memset instead of a for loop

   memset(marked,0 , CacheBlock_size);

   This had a substantial effect, reducing runtime from $12.7$s → $9.7$s

   

2. First attempt to Unroll Loops

   Loop control introduces conditional branches. Processing several elements per iteration can reduce the overhead and branch prediction work.

   for (i = first; i < CacheBlock_size; i += 4*prime) {
       marked[i] = 1;
       if (i + prime < CacheBlock_size) marked[i + prime] = 1;
       if (i + 2*prime < CacheBlock_size) marked[i + 2*prime] = 1;
       if (i + 3*prime < CacheBlock_size) marked[i + 3*prime] = 1;
   }
   

   This version marks four marked entries per iteration, bringing runtime to $7.7$ s.

   

   The counting loop can be unrolled in the same way:

       for (i = 0; i < CacheBlock_size; i += 4) {
           if (!marked[i]) count_cacheBlock++;
           if (i + 1 < CacheBlock_size && !marked[i + 1]) count_cacheBlock++;
           if (i + 2 < CacheBlock_size && !marked[i + 2]) count_cacheBlock++;
           if (i + 3 < CacheBlock_size && !marked[i + 3]) count_cacheBlock++;
       }
   

   Runtime dropped to $6.1$ s.

   The loop above still contains many if statements. Since the point of unrolling is to reduce branch overhead, I adjusted the structure and the values stored in marked:

           for (i = 0; i < CacheBlock_size; i += 4) {
               count_cacheBlock0 += marked[i];
               if (i + 3 < CacheBlock_size) {
                   count_cacheBlock1 += marked[i + 1];
                   count_cacheBlock2 += marked[i + 2];
                   count_cacheBlock3 += marked[i + 3];
               } else if (i + 2 < CacheBlock_size) {
                   count_cacheBlock1 += marked[i + 1];
                   count_cacheBlock2 += marked[i + 2];
               } else {
                   count_cacheBlock1 += marked[i + 1];
               }
           }
   

   Now marked is initialized to 1, and composite entries are set to 0, so the entries can be summed directly. Runtime dropped to $5.3$ s.

   

   Actually we can remove more if branches. Will be introduced later.

3. Count with __builtin_popcountll

   __builtin_popcountll is a GCC built-in function that counts the 1 bits in an unsigned long long value, also known as its Hamming weight.

   tmpSize = (CacheBlock_size >> 6) << 6;
   for (i = 0; i < tmpSize; i += 64) {
       count_cacheBlock0 += __builtin_popcountll(*((unsigned long long *)(marked + i)));
       count_cacheBlock1 += __builtin_popcountll(*((unsigned long long *)(marked + i + 8)));
       count_cacheBlock2 += __builtin_popcountll(*((unsigned long long *)(marked + i + 16)));
       count_cacheBlock3 += __builtin_popcountll(*((unsigned long long *)(marked + i + 24)));
       count_cacheBlock4 += __builtin_popcountll(*((unsigned long long *)(marked + i + 32)));
       count_cacheBlock5 += __builtin_popcountll(*((unsigned long long *)(marked + i + 40)));
       count_cacheBlock6 += __builtin_popcountll(*((unsigned long long *)(marked + i + 48)));
       count_cacheBlock7 += __builtin_popcountll(*((unsigned long long *)(marked + i + 56)));
   }
   count += count_cacheBlock0 + count_cacheBlock1 + count_cacheBlock2 + count_cacheBlock3 + count_cacheBlock4 + count_cacheBlock5 + count_cacheBlock6 + count_cacheBlock7;
   for (; i < CacheBlock_size; i++)
       count += marked[i];
   

   Each __builtin_popcountll call counts the set bits across eight byte-sized marked entries. With eight calls, an iteration handles 64 entries, which substantially improves efficiency. To avoid reading past the array boundary, first clear the lowest six bits of CacheBlock_size and handle the remaining entries separately.

   Runtime reached $4.3$ s. The reported prime count was incorrect because I forgot to reset the counters; the source reports that this did not affect the timing.

4. Further Reduce Branches in the Marking Loop

   Calculate a safe limit for the unrolled loop, then handle the remaining entries separately:

   tmpSize = CacheBlock_size - (CacheBlock_size - first) % (prime << 2);
   for (i = first; i < tmpSize; i += (prime << 2)) {
       marked[i] = 0;
       marked[i + prime] = 0;
       marked[i + (prime << 1)] = 0;
       marked[i + 3 * prime] = 0;
   }
   for (; i < CacheBlock_size; i += prime) marked[i] = 0;
   

   Runtime dropped to $3.5$ s.

5. Increase the Marking Loop's Unroll Factor

   Instead of handling four entries per iteration, handle eight:

   tmpSize = CacheBlock_size - (CacheBlock_size - first) % (prime << 3);
   for (i = first; i < tmpSize; i += (prime << 3)) {
       marked[i] = 0;
       marked[i + prime] = 0;
       marked[i + (prime << 1)] = 0;
       marked[i + 3 * prime] = 0;
       marked[i + (prime << 2)] = 0;
       marked[i + 5 * prime] = 0;
       marked[i + 6 * prime] = 0;
       marked[i + 7 * prime] = 0;
   }
   for (; i < CacheBlock_size; i += prime) marked[i] = 0;
   

   This was slightly faster, reaching $3.3$ s. A change this small might simply reflect normal variation, but at least it did not make performance worse.

6. Miscellaneous Optimizations

   • Use pointer operations:

     char *mark = marked + i;
     *mark = 0;
     *(mark + prime) = 0;
     *(mark + (prime << 1)) = 0;
     ...
     

   • Use more bitwise operations, such as replacing %2 with &1. The source describes bitwise operations as much faster than ordinary arithmetic.

   • Use __builtin_expect to provide branch prediction hints. For branches with several alternatives, I counted how often each path was taken and added likely(x) or unlikely(x) accordingly. This took considerable effort, but the improvement was modest.

     #define likely(x) __builtin_expect(!!(x), 1)
     #define unlikely(x) __builtin_expect(!!(x), 0)
     

   • Store Block_low_value % prime in a variable. This expression occurs repeatedly in the loop body, and avoiding repeated calculation improved performance.

   • Adjust CacheBlock_size, as discussed under optimize3. This produced a substantial improvement.

7. Use register and Inline Assembly

   The register keyword suggests storing a variable in a CPU register for faster access. I initially expected little benefit because the keyword is obsolete and modern compilers may ignore the hint. Surprisingly, it helped in these tests. After several experiments, I added register to all variables used in the loop body. With 8 processes and a problem size of $4\times10^9$, this reduced runtime by about 1 second—a substantial improvement after many other experiments had yielded little progress.

   I found that register was deprecated in C++11 and removed in C++17, so I checked the platform's mpic++ version:

   

   The GCC 7.5.0 documentation indicates that G++ 7.5 uses the GNU dialect of C++14 by default.

   GNU dialect of -std=c++14. This is the default for C++ code.

   The C++ dialect options documentation explains that the compiler can warn about the register storage class specifier, except when used with GNU explicit register variables. It also records its deprecation in C++11 and removal in C++17, with the warning enabled by default under -std=c++1z.

   Warn on uses of the register storage class specifier, except when it is part of the GNU Explicit Register Variables extension. The use of the register keyword as storage class specifier has been deprecated in C++11 and removed in C++17. Enabled by default with -std=c++1z.

   The variables in this code appear to be local variables, so I continued to Local Register Variables. This describes explicit register assignments using asm and their use with inline assembly:

   register int *p1 asm ("r0") = …;
   register int *p2 asm ("r1") = …;
   register int *result asm ("r0");
   asm ("sysint" : "=r" (result) : "0" (p1), "r" (p2));
   

   Following that pattern, I rewrote the marking loop in inline assembly while retaining loop unrolling:

       register bool *tmp_marked asm("x7");
       register unsigned LL tmp_i asm("x9");
       register unsigned LL tmp_prime asm("x10");
       register unsigned int zero asm("w8");
   
               tmpSize = CacheBlock_size - (CacheBlock_size - first) % (prime << 3);
               i = first;
               tmp_marked = &marked[0];
               tmp_i = (unsigned LL)i;
               tmp_prime = (unsigned LL)prime;
               zero = 0;
               asm volatile(
                   "2:\n"
                   "cmp x9, %2\n"       // i < tmpSize?
                   "bge 1f\n"
                   "strb w8, [x7,x9]\n" // marked[i] = 0
                   "adds x9, x9, x10\n" // i+=prime
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "strb w8, [x7,x9]\n"
                   "adds x9, x9, x10\n"
                   "b 2b\n"
                   "1:\n"
                   : "+r"(tmp_i)
                   : "r"(tmp_marked), "r"(tmpSize), "r"(zero)
                   : "cc", "memory");
               i = tmp_i;
               for (; i < CacheBlock_size; i += prime) marked[i] = 0;
   

   Runtime improved slightly to $2.7$ s.

   I also tried replacing the code that counts 1s in marked with inline assembly, but the results were disappointing. Perhaps __builtin_popcount was already fast enough, or perhaps my assembly implementation was insufficient.

Results for optimize4.cpp:

With 1–16 processes and a problem size of $3\times10^9$:

Final evaluation configuration: 8 processes and a problem size of $4\times10^9$.

Source Code

base.cpp

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define LL long long

int main(int argc, char *argv[]) {
    LL count;            /* Local prime count */
    double elapsed_time; /* Parallel execution time */
    LL first;            /* Index of first multiple */
    LL global_count = 0; /* Global prime count */
    LL high_value;       /* Highest value on this proc */
    LL i;
    int id;              /* Process ID number */
    LL index;            /* Index of current prime */
    LL low_value;        /* Lowest value on this proc */
    char *marked;        /* Portion of 2,...,'n' */
    LL n;                /* Sieving from 2, ..., 'n' */
    int p;               /* Number of processes */
    LL proc0_size;       /* Size of proc 0's subarray */
    LL prime;            /* Current prime */
    LL size;             /* Elements in 'marked' */

    MPI_Init(&argc, &argv);

    /* Start the timer */

    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id) printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }

    n = atoll(argv[1]);

    /* Figure out this process's share of the array, as
        well as the integers represented by the first and
       last array elements */

    low_value = 2 + id * (n - 1) / p;
    high_value = 1 + (id + 1) * (n - 1) / p;
    size = high_value - low_value + 1;

    /* Bail out if all the primes used for sieving are
       not all held by process 0 */

    proc0_size = (n - 1) / p;

    if ((2 + proc0_size) < (int)sqrt((double)n)) {
        if (!id) printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }

    /* Allocate this process's share of the array. */

    marked = (char *)malloc(size);

    if (marked == NULL) {
        printf("Cannot allocate enough memory\n");
        MPI_Finalize();
        exit(1);
    }

    //!------这里也许可以用 memset
    for (i = 0; i < size; i++) marked[i] = 0;
    if (!id) index = 0;
    prime = 2;
    do {
        if (prime * prime > low_value)
            first = prime * prime - low_value;
        else {
            if (!(low_value % prime))
                first = 0;
            else
                first = prime - (low_value % prime);
        }
        for (i = first; i < size; i += prime) marked[i] = 1;
        if (!id) {
            while (marked[++index])
                ;
            prime = index + 2;
        }
        if (p > 1) MPI_Bcast(&prime, 1, MPI_INT, 0, MPI_COMM_WORLD);
    } while (prime * prime <= n);
    count = 0;
    for (i = 0; i < size; i++)
        if (!marked[i]) count++;
    MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM,
               0, MPI_COMM_WORLD);

    /* Stop the timer */

    elapsed_time += MPI_Wtime();

    /* Print the results */

    if (!id) {
        printf("There are %lld primes less than or equal to %lld\n",
               global_count, n);
        printf("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize();
    return 0;
}

optimize.cpp

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define LL long long

int main(int argc, char *argv[]) {
    LL count;            /* Local prime count */
    double elapsed_time; /* Parallel execution time */
    LL first;            /* Index of first multiple */
    LL global_count = 0; /* Global prime count */
    LL high_value;       /* Highest value on this proc */
    LL i;
    int id;              /* Process ID number */
    LL index;            /* Index of current prime */
    LL low_value;        /* Lowest value on this proc */
    char *marked;        /* Portion of 2,...,'n' */
    LL n;                /* Sieving from 2, ..., 'n' */
    int p;               /* Number of processes */
    LL proc0_size;       /* Size of proc 0's subarray */
    LL prime;            /* Current prime */
    LL size;             /* Elements in 'marked' */

    MPI_Init(&argc, &argv);

    /* Start the timer */

    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id)
            printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }

    n = atoll(argv[1]);

    /* Figure out this process's share of the array, as
       well as the integers represented by the first and
       last array elements */

    low_value = 2 + id * (n - 1) / p;
    high_value = 1 + (id + 1) * (n - 1) / p;
    // 大小开一半就行
    size = (high_value - low_value) / 2 + 1;

    /* Bail out if all the primes used for sieving are
       not all held by process 0 */

    proc0_size = (n - 1) / p;

    if ((2 + proc0_size) < (int)sqrt((double)n)) {
        if (!id)
            printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }

    /* Allocate this process's share of the array. */

    marked = (char *)malloc(size);

    if (marked == NULL) {
        printf("Cannot allocate enough memory\n");
        MPI_Finalize();
        exit(1);
    }

    for (i = 0; i < size; i++)
        marked[i] = 0;
    if (!id)
        index = 0;
    prime = 3;
    do {
        if (prime * prime > low_value)
            first = (prime * prime - low_value) / 2;
        else {
            if (!(low_value % prime))
                first = 0;
            else if (low_value % prime % 2 == 0)
                first = prime - ((low_value % prime) / 2);
            else
                first = (prime - (low_value % prime)) / 2;
        }
        for (i = first; i < size; i += prime)
            marked[i] = 1;
        if (!id) {
            while (marked[++index])
                ;
            prime = index * 2 + 3;
        }
        if (p > 1)
            MPI_Bcast(&prime, 1, MPI_INT, 0, MPI_COMM_WORLD);
    } while (prime * prime <= n);
    count = 0;
    for (i = 0; i < size; i++)
        if (!marked[i])
            count++;

    MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
    global_count++;
    /* Stop the timer */

    elapsed_time += MPI_Wtime();

    /* Print the results */

    if (!id) {
        printf("There are %lld primes less than or equal to %lld\n", global_count,
               n);
        printf("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize();
    return 0;
}

optimize2.cpp

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
using namespace std;
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define LL long long

int main(int argc, char *argv[]) {
    LL count;            /* Local prime count */
    double elapsed_time; /* Parallel execution time */
    LL first;            /* Index of first multiple */
    LL global_count = 0; /* Global prime count */
    LL high_value;       /* Highest value on this proc */
    LL i;
    int id;              /* Process ID number */
    LL index;            /* Index of current prime */
    LL low_value;        /* Lowest value on this proc */
    char *marked;        /* Portion of 2,...,'n' */
    LL n;                /* Sieving from 2, ..., 'n' */
    int p;               /* Number of processes */
    LL proc0_size;       /* Size of proc 0's subarray */
    LL prime;            /* Current prime */
    LL size;             /* Elements in 'marked' */

    MPI_Init(&argc, &argv);

    /* Start the timer */

    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id)
            printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }

    n = atoll(argv[1]);

    // 欧拉筛sqrt(n)内的质数
    vector<int> pri;
    int pri_i = 1;
    int sqrt_N = sqrt(n) + 10;
    bool not_prime[sqrt_N + 10] = {0};
    for (int i = 2; i <= sqrt_N; i++) {
        if (!not_prime[i]) pri.push_back(i);
        for (int prime_i : pri) {
            if (i * prime_i > sqrt_N) break;
            not_prime[i * prime_i] = true;
            if (i % prime_i == 0) break;
        }
    }

    /* Figure out this process's share of the array, as
       well as the integers represented by the first and
       last array elements */

    low_value = 2 + id * (n - 1) / p;
    high_value = 1 + (id + 1) * (n - 1) / p;
    // 大小开一半就行
    size = (high_value - low_value) / 2 + 1;

    /* Bail out if all the primes used for sieving are
       not all held by process 0 */

    proc0_size = (n - 1) / p;

    if ((2 + proc0_size) < (int)sqrt((double)n)) {
        if (!id)
            printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }

    /* Allocate this process's share of the array. */

    marked = (char *)malloc(size);

    if (marked == NULL) {
        printf("Cannot allocate enough memory\n");
        MPI_Finalize();
        exit(1);
    }

    for (i = 0; i < size; i++)
        marked[i] = 0;
    if (!id)
        index = 0;
    prime = pri[pri_i];
    do {
        if (prime * prime > low_value)
            first = (prime * prime - low_value) / 2;
        else {
            if (!(low_value % prime))
                first = 0;
            else if (low_value % prime % 2 == 0) // 去除偶数
                first = prime - ((low_value % prime) / 2);
            else
                first = (prime - (low_value % prime)) / 2;
        }
        for (i = first; i < size; i += prime)
            marked[i] = 1;
        prime = pri[++pri_i];
    } while (prime * prime <= n);
    count = 0;
    for (i = 0; i < size; i++)
        if (!marked[i])
            count++;

    MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
    global_count++;
    /* Stop the timer */

    elapsed_time += MPI_Wtime();

    /* Print the results */

    if (!id) {
        printf("There are %lld primes less than or equal to %lld\n", global_count,
               n);
        printf("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize();
    return 0;
}

optimize3.cpp

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
using namespace std;
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define LL long long
#define CACHE_SIZE 33554432
#define CACHELINE_SIZE 64

int main(int argc, char *argv[]) {
    LL count;            /* Local prime count */
    double elapsed_time; /* Parallel execution time */
    LL first;            /* Index of first multiple */
    LL global_count = 0; /* Global prime count */
    LL high_value;       /* Highest value on this proc */
    LL i;
    int id;              /* Process ID number */
    LL index;            /* Index of current prime */
    LL low_value;        /* Lowest value on this proc */
    char *marked;        /* Portion of 2,...,'n' */
    LL n;                /* Sieving from 2, ..., 'n' */
    int p;               /* Number of processes */
    LL proc0_size;       /* Size of proc 0's subarray */
    LL prime;            /* Current prime */
    LL size;             /* Elements in 'marked' */

    MPI_Init(&argc, &argv);

    /* Start the timer */

    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id)
            printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }

    n = atoll(argv[1]);

    // 欧拉筛sqrt(n)内的质数
    vector<int> pri;
    int pri_i = 1;
    int sqrt_N = sqrt(n) + 10;
    bool not_prime[sqrt_N + 10] = {0};
    for (int i = 2; i <= sqrt_N; i++) {
        if (!not_prime[i]) pri.push_back(i);
        for (int prime_i : pri) {
            if (i * prime_i > sqrt_N) break;
            not_prime[i * prime_i] = true;
            if (i % prime_i == 0) break;
        }
    }

    /* Figure out this process's share of the array, as
       well as the integers represented by the first and
       last array elements */
    LL N = (n - 1) / 2;
    LL low_index = id * (N / p) + MIN(id, N % p); // 该进程的第一个数的索引
    LL high_index = (id + 1) * (N / p) + MIN(id + 1, N % p) - 1;
    low_value = low_index * 2 + 3;
    high_value = (high_index + 1) * 2 + 1;

    size = (high_value - low_value) / 2 + 1;

    /* Bail out if all the primes used for sieving are
       not all held by process 0 */

    proc0_size = (n - 1) / p;

    if ((2 + proc0_size) < (int)sqrt((double)n)) {
        if (!id)
            printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }

    /* Allocate this process's share of the array. */
    //

    int CacheBlock_size = (CACHE_SIZE / p) >> 4;
    int Block_need = size / CacheBlock_size;
    int Block_rest = size % CacheBlock_size;
    int Block_id = 0;
    int Block_N = CacheBlock_size - 1;
    LL Block_low_value, Block_high_value, Block_first, Block_last;
    LL count_cacheBlock;

    marked = (char *)malloc(CacheBlock_size);

    if (marked == NULL) {
        printf("Cannot allocate enough memory\n");
        MPI_Finalize();
        exit(1);
    }

    count = 0;
    while (Block_id <= Block_need) {
        // 更新Cache
        Block_first = Block_id * Block_N + MIN(Block_id, Block_rest) + low_index;
        Block_last = (Block_id + 1) * Block_N + MIN(Block_id + 1, Block_rest) - 1 + low_index;
        Block_low_value = 2 * Block_first + 3;
        if (Block_id == Block_need) {
            Block_high_value = high_value;
            Block_last = high_index;
            CacheBlock_size = (Block_high_value - Block_low_value) / 2 + 1;
        } else {
            Block_high_value = 2 * Block_last + 3;
        }

        pri_i = 1;
        prime = pri[pri_i];
        count_cacheBlock = 0;
        for (int i = 0; i < CacheBlock_size; i++) marked[i] = 0;

        do {
            if (prime * prime > Block_low_value)
                first = (prime * prime - Block_low_value) / 2;
            else {
                if (!(Block_low_value % prime))
                    first = 0;
                else if (Block_low_value % prime % 2 == 0)
                    first = prime - (Block_low_value % prime) / 2;
                else
                    first = (prime - (Block_low_value % prime)) / 2;
            }
            for (i = first; i < CacheBlock_size; i += prime)
                marked[i] = 1;

            prime = pri[++pri_i];
        } while (prime * prime <= Block_high_value);

        for (i = 0; i < CacheBlock_size; i++)
            if (!marked[i])
                count_cacheBlock++;

        count += count_cacheBlock;
        Block_id++;
    }

    MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
    global_count++;
    /* Stop the timer */

    elapsed_time += MPI_Wtime();

    /* Print the results */

    if (!id) {
        printf("There are %lld primes less than or equal to %lld\n", global_count,
               n);
        printf("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize();
    return 0;
}

optimize4.cpp

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <cstring>
using namespace std;
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define LL long long
#define CACHE_SIZE 33554432
#define likely(x) __builtin_expect(!!(x), 1)
#define unlikely(x) __builtin_expect(!!(x), 0)

int main(int argc, char *argv[]) {
    LL count = 0;           /* Local prime count */
    double elapsed_time;    /* Parallel execution time */
    register LL first;      /* Index of first multiple */
    LL global_count = 0;    /* Global prime count */
    register LL high_value; /* Highest value on this proc */
    register LL i;
    register int id;        /* Process ID number */
    LL index;               /* Index of current prime */
    LL low_value;           /* Lowest value on this proc */
    register bool *marked;  /* Portion of 2,...,'n' */
    register bool *mark;
    LL n;                   /* Sieving from 2, ..., 'n' */
    int p;                  /* Number of processes */
    LL proc0_size;          /* Size of proc 0's subarray */
    register LL prime;      /* Current prime */
    LL size;                /* Elements in 'marked' */
    register LL Block_low_value, Block_high_value, Block_first, Block_last;
    register LL prime_square;
    register LL tmpSize; // 用来循环计数的，防止越界
    register LL tmpPrime;
    register LL count_cacheBlock0 = 0, count_cacheBlock1 = 0, count_cacheBlock2 = 0, count_cacheBlock3 = 0, count_cacheBlock4 = 0, count_cacheBlock5 = 0, count_cacheBlock6 = 0, count_cacheBlock7 = 0;

    register bool *tmp_marked asm("x7");
    register unsigned LL tmp_i asm("x9");
    register unsigned LL tmp_prime asm("x10");
    register unsigned int zero asm("w8");

    MPI_Init(&argc, &argv);

    /* Start the timer */

    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id)
            printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }

    n = atoll(argv[1]);

    // 欧拉筛sqrt(n)内的质数

    register vector<int> pri;
    int pri_i = 1;
    register int sqrt_N = sqrt(n) + 10;
    register bool not_prime[sqrt_N + 10] = {0};
    for (register int i = 2; i <= sqrt_N; i++) {
        if (!not_prime[i]) pri.emplace_back(i);
        for (register int prime_i : pri) {
            if (i * prime_i > sqrt_N) break;
            not_prime[i * prime_i] = true;
            if (i % prime_i == 0) break;
        }
    }

    /* Figure out this process's share of the array, as
       well as the integers represented by the first and
       last array elements */
    LL N = (n - 1) >> 1;
    register LL low_index = id * (N / p) + MIN(id, N % p); // 该进程的第一个数的索引
    LL high_index = (id + 1) * (N / p) + MIN(id + 1, N % p) - 1;
    low_value = (low_index << 1) + 3;
    high_value = ((high_index + 1) << 1) + 1;
    size = ((high_value - low_value) >> 1) + 1;

    /* Bail out if all the primes used for sieving are
       not all held by process 0 */

    proc0_size = (n - 1) / p;

    if ((2 + proc0_size) < (int)sqrt((double)n)) {
        if (!id)
            printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }

    /* Allocate this process's share of the array. */
    //

    int CacheBlock_size = (CACHE_SIZE / p) >> 4;
    register int Block_need = size / CacheBlock_size;
    register int Block_rest = size % CacheBlock_size;
    register int Block_id = 0;
    register int Block_N = CacheBlock_size - 1;
    register int BlocLowModPrime;
    marked = (bool *)malloc(CacheBlock_size);

    if (marked == NULL) {
        printf("Cannot allocate enough memory\n");
        MPI_Finalize();
        exit(1);
    }

    while (Block_id <= Block_need) {
        pri_i = 1;
        prime = 3;
        prime_square = prime * prime;
        memset(marked, 1, CacheBlock_size);
        count_cacheBlock0 = 0, count_cacheBlock1 = 0, count_cacheBlock2 = 0, count_cacheBlock3 = 0, count_cacheBlock4 = 0, count_cacheBlock5 = 0, count_cacheBlock6 = 0, count_cacheBlock7 = 0;
        // 更新Cache
        Block_first = Block_id * Block_N + MIN(Block_id, Block_rest) + low_index;
        Block_last = (Block_id + 1) * Block_N + MIN(Block_id + 1, Block_rest) - 1 + low_index;
        Block_low_value = (Block_first << 1) + 3;
        if (unlikely(Block_id == Block_need)) {
            Block_high_value = high_value;
            Block_last = high_index;
            CacheBlock_size = ((Block_high_value - Block_low_value) >> 1) + 1;
        } else {
            Block_high_value = (Block_last << 1) + 3;
        }

        do {
            BlocLowModPrime = Block_low_value % prime;
            if (unlikely(prime_square > Block_low_value))
                first = (prime_square - Block_low_value) >> 1;
            else {
                if (unlikely(BlocLowModPrime == 0))
                    first = 0;
                else {
                    first = ((BlocLowModPrime & 1) == 0) ?
                                (prime - (BlocLowModPrime >> 1)) :
                                (prime - BlocLowModPrime) >> 1;
                }
            }

            tmpSize = CacheBlock_size - (CacheBlock_size - first) % (prime << 3);
            i = first;
            tmp_marked = &marked[0];        // x7
            tmp_i = (unsigned LL)i;         // x9
            tmp_prime = (unsigned LL)prime; // x10
            zero = 0;                       // w8
            asm volatile(
                "2:\n"
                "cmp x9, %2\n"       // i < tmpSize?
                "bge 1f\n"
                "strb w8, [x7,x9]\n" // marked[i] = 0
                "adds x9, x9, x10\n" // i+=prime
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "strb w8, [x7,x9]\n"
                "adds x9, x9, x10\n"
                "b 2b\n"
                "1:\n"
                : "+r"(tmp_i)
                : "r"(tmp_marked), "r"(tmpSize), "r"(zero)
                : "cc", "memory");
            i = tmp_i;
            for (; i < CacheBlock_size; i += prime) marked[i] = 0;

            prime = pri[++pri_i];
            prime_square = prime * prime;
        } while (prime_square <= Block_high_value);

        tmpSize = (CacheBlock_size >> 6) << 6;
        for (i = 0; i < tmpSize; i += 64) {
            mark = marked + i;
            count_cacheBlock0 += __builtin_popcountll(*((unsigned long long *)(mark)));
            count_cacheBlock1 += __builtin_popcountll(*((unsigned long long *)(mark + 8)));
            count_cacheBlock2 += __builtin_popcountll(*((unsigned long long *)(mark + 16)));
            count_cacheBlock3 += __builtin_popcountll(*((unsigned long long *)(mark + 24)));
            count_cacheBlock4 += __builtin_popcountll(*((unsigned long long *)(mark + 32)));
            count_cacheBlock5 += __builtin_popcountll(*((unsigned long long *)(mark + 40)));
            count_cacheBlock6 += __builtin_popcountll(*((unsigned long long *)(mark + 48)));
            count_cacheBlock7 += __builtin_popcountll(*((unsigned long long *)(mark + 56)));
        }
        count += count_cacheBlock0 + count_cacheBlock1 + count_cacheBlock2 + count_cacheBlock3 + count_cacheBlock4 + count_cacheBlock5 + count_cacheBlock6 + count_cacheBlock7;
        for (; i < CacheBlock_size; i++) count += marked[i];

        Block_id++;
    }

    MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
    global_count++;

    /* Stop the timer */

    elapsed_time += MPI_Wtime();

    /* Print the results */

    if (!id) {
        printf("There are %lld primes less than or equal to %lld\n", global_count,
               n);
        printf("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize();
    return 0;
}