Parallel Algorithm Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Parallel exclusive scan using two-phase block method.
5
// Phase 1: Each thread computes exclusive scan on its block and records block sum.
6
// Phase 2: Compute prefix over block sums, then add offset to each block in parallel.
7
8
vector<long long> parallel_exclusive_scan(const vector<int>& a, unsigned num_threads = thread::hardware_concurrency()) {
9
    const size_t n = a.size();
10
    vector<long long> out(n, 0);
11
    if (n == 0) return out;
12
13
    if (num_threads == 0) num_threads = 4; // fallback
14
    num_threads = min<unsigned>(num_threads, static_cast<unsigned>(n));
15
    const size_t block_size = (n + num_threads - 1) / num_threads;
16
17
    vector<long long> blockSums(num_threads, 0);
18
    vector<thread> threads;
19
    threads.reserve(num_threads);
20
21
    // Phase 1: local exclusive scans per block
22
    for (unsigned t = 0; t < num_threads; ++t) {
23
        threads.emplace_back([&, t](){
24
            const size_t begin = t * block_size;
25
            const size_t end = min(n, begin + block_size);
26
            long long sum = 0;
27
            for (size_t i = begin; i < end; ++i) {
28
                out[i] = sum;       // exclusive: write previous sum
29
                sum += a[i];        // accumulate current element
30
            }
31
            blockSums[t] = sum;     // total of this block
32
        });
33
    }
34
    for (auto &th : threads) th.join();
35
36
    // Compute offsets as exclusive prefix of blockSums (sequential, size = num_threads)
37
    vector<long long> offsets(num_threads, 0);
38
    long long acc = 0;
39
    for (unsigned t = 0; t < num_threads; ++t) {
40
        offsets[t] = acc;
41
        acc += blockSums[t];
42
    }
43
44
    // Phase 2: add offsets in parallel
45
    threads.clear();
46
    for (unsigned t = 0; t < num_threads; ++t) {
47
        threads.emplace_back([&, t](){
48
            const size_t begin = t * block_size;
49
            const size_t end = min(n, begin + block_size);
50
            const long long off = offsets[t];
51
            for (size_t i = begin; i < end; ++i) {
52
                out[i] += off;
53
            }
54
        });
55
    }
56
    for (auto &th : threads) th.join();
57
58
    return out;
59
}
60
61
int main() {
62
    ios::sync_with_stdio(false);
63
    cin.tie(nullptr);
64
65
    vector<int> a = {3, 1, 7, 0, 4, 1, 6, 3};
66
    auto res = parallel_exclusive_scan(a);
67
68
    // Verify against sequential exclusive scan
69
    vector<long long> check(a.size(), 0);
70
    long long s = 0;
71
    for (size_t i = 0; i < a.size(); ++i) {
72
        check[i] = s;
73
        s += a[i];
74
    }
75
76
    cout << "Input:       "; for (auto v: a) cout << v << ' '; cout << '\n';
77
    cout << "Parallel out:"; for (auto v: res) cout << ' ' << v; cout << '\n';
78
    cout << "Sequential :"; for (auto v: check) cout << ' ' << v; cout << '\n';
79
80
    bool ok = (res == check);
81
    cout << (ok ? "Match" : "Mismatch") << '\n';
82
    return ok ? 0 : 1;
83
}
84

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Parallel merge sort using std::async with a cutoff to limit task overhead.
5
// The algorithm remains work-efficient (O(n log n)) while reducing span through concurrency.
6
7
static const size_t CUTOFF = 1 << 14; // switch to sequential sort for small arrays
8
9
void merge_ranges(vector<int>& a, vector<int>& tmp, size_t l, size_t m, size_t r) {
10
    size_t i = l, j = m, k = l;
11
    while (i < m && j < r) tmp[k++] = (a[i] <= a[j]) ? a[i++] : a[j++];
12
    while (i < m) tmp[k++] = a[i++];
13
    while (j < r) tmp[k++] = a[j++];
14
    // copy back
15
    for (size_t t = l; t < r; ++t) a[t] = tmp[t];
16
}
17
18
void parallel_merge_sort_rec(vector<int>& a, vector<int>& tmp, size_t l, size_t r, unsigned depth) {
19
    const size_t n = r - l;
20
    if (n <= 1) return;
21
    if (n <= CUTOFF || depth == 0) {
22
        sort(a.begin() + l, a.begin() + r);
23
        return;
24
    }
25
    size_t m = l + n / 2;
26
    // Spawn one half asynchronously and sort the other in current thread
27
    auto fut = async(launch::async, [&]{ parallel_merge_sort_rec(a, tmp, l, m, depth - 1); });
28
    parallel_merge_sort_rec(a, tmp, m, r, depth - 1);
29
    fut.get();
30
    merge_ranges(a, tmp, l, m, r);
31
}
32
33
void parallel_merge_sort(vector<int>& a) {
34
    vector<int> tmp(a.size());
35
    unsigned hw = thread::hardware_concurrency();
36
    if (hw == 0) hw = 4;
37
    // Depth controls number of concurrent tasks: up to 2^depth tasks
38
    unsigned depth = static_cast<unsigned>(log2(hw));
39
    parallel_merge_sort_rec(a, tmp, 0, a.size(), depth);
40
}
41
42
int main() {
43
    ios::sync_with_stdio(false);
44
    cin.tie(nullptr);
45
46
    vector<int> a(1 << 20);
47
    iota(a.begin(), a.end(), 0);
48
    // Shuffle to create a random permutation
49
    mt19937 rng(123);
50
    shuffle(a.begin(), a.end(), rng);
51
52
    vector<int> b = a;
53
    parallel_merge_sort(a);
54
    sort(b.begin(), b.end());
55
56
    cout << (a == b ? "Sorted correctly" : "Incorrectly sorted") << "\n";
57
    return a == b ? 0 : 1;
58
}
59

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute C = A * B where A, B, C are n x n matrices (dense), parallelized by rows.
5
// For better cache locality, we transpose B first so that inner loop reads contiguous memory.
6
7
typedef vector<double> Vec;
8
typedef vector<Vec> Mat;
9
10
Mat transpose(const Mat& B) {
11
    size_t n = B.size();
12
    Mat Bt(n, Vec(n));
13
    for (size_t i = 0; i < n; ++i)
14
        for (size_t j = 0; j < n; ++j)
15
            Bt[j][i] = B[i][j];
16
    return Bt;
17
}
18
19
Mat parallel_matmul(const Mat& A, const Mat& B, unsigned num_threads = thread::hardware_concurrency()) {
20
    size_t n = A.size();
21
    Mat C(n, Vec(n, 0.0));
22
    if (n == 0) return C;
23
    if (num_threads == 0) num_threads = 4;
24
    num_threads = min<unsigned>(num_threads, static_cast<unsigned>(n));
25
26
    Mat Bt = transpose(B); // improves memory access patterns
27
28
    size_t rows_per_thread = (n + num_threads - 1) / num_threads;
29
    vector<thread> threads;
30
    threads.reserve(num_threads);
31
32
    for (unsigned t = 0; t < num_threads; ++t) {
33
        threads.emplace_back([&, t]() {
34
            size_t r_begin = t * rows_per_thread;
35
            size_t r_end = min(n, r_begin + rows_per_thread);
36
            for (size_t i = r_begin; i < r_end; ++i) {
37
                for (size_t j = 0; j < n; ++j) {
38
                    const double* ai = A[i].data();
39
                    const double* btj = Bt[j].data();
40
                    double sum = 0.0;
41
                    // Inner product of A[i][:] and B[:][j] (via Bt[j][:])
42
                    for (size_t k = 0; k < n; ++k) sum += ai[k] * btj[k];
43
                    C[i][j] = sum;
44
                }
45
            }
46
        });
47
    }
48
    for (auto& th : threads) th.join();
49
50
    return C;
51
}
52
53
int main() {
54
    ios::sync_with_stdio(false);
55
    cin.tie(nullptr);
56
57
    size_t n = 4;
58
    Mat A(n, Vec(n)), B(n, Vec(n));
59
    // Simple test: A = identity, C should equal B
60
    for (size_t i = 0; i < n; ++i) {
61
        for (size_t j = 0; j < n; ++j) {
62
            A[i][j] = (i == j) ? 1.0 : 0.0;
63
            B[i][j] = (i + j);
64
        }
65
    }
66
    Mat C = parallel_matmul(A, B);
67
    bool ok = true;
68
    for (size_t i = 0; i < n; ++i)
69
        for (size_t j = 0; j < n; ++j)
70
            ok = ok && (fabs(C[i][j] - B[i][j]) < 1e-9);
71
72
    cout << (ok ? "Matmul correct" : "Matmul incorrect") << "\n";
73
    return ok ? 0 : 1;
74
}
75