^-^
Assignment for Parallel Computing
  1. assignment 01
    1. 1 Work and Span
    2. 2 Y++ Language
    3. 3 Compute the Integral
  2. Experiment 1: Parallel Eratosthenes Sieve (MPI)
    1. MPI Something
    2. 并行思路
    3. 优化
    4. 更多优化思路
  3. Parallel Floyd Algorithm
    1. MPI 点对点通信

assignment 01

1 Work and Span

  1. 一个二重循环,

  2. 递推关系: . 因此 .

  3. . 解得 .

  4. , , 则 .

  6. , , 则 ,

  1. , ,

  1. , .

2 Y++ Language

给定一个仅含有 的括号序列判断是否匹配.

  1. 设计一个时间复杂度为 的串行算法, 要求只能使用 的额外的空间.

SOL 首先想到的做法是维护一个栈, 从括号序列的结尾开始遍历, 如果是右括号则入栈; 如果是左括号则弹出栈顶元素, 得到一对匹配的括号, 当串空时栈也空, 该括号序列匹配, 反之不匹配. 这个要求额外开 的栈空间, 不满足题意.

但是借此思路可以不存括号而是对括号计数, 即定义右括号数目的变量, 从括号序列的结尾开始遍历, 遇到右括号自增, 遇到左括号自减, 串空时计数变量为 0 则匹配, 反之不匹配.

  1. 设计用来计数一个括号序列中左括号 ( 和 右括号 ) 的个数的并行算法, 可以使用额外的 的空间, 保证 的 work 和 的 span.

SOL 类似于 reduce, 将整个序列分成两部分分别并行地计数其左右括号数, 再 join 起来即可.

1
2
3
4
5
6
7
8
9
10
11
12
def count(P,n):
    if n==0:
        return [0,0]
    if n==1:
        if P[0]='(':
            return [1,0]
        else:
            return [0,1]
    par_do:
        cnt1=count(P,n/2)
        cnt2=count(P+n/2,n-n/2)
   	return [cnt1[0]+cnt2[0],cnt1[1]+cnt2[1]]

满足 , .

  1. 将 (1) 中的算法并行化, 可以使用额外的 的空间, 保证 的 work 和 的 span.

SOL 分支, 将序列分成两份, 分别用一个计数变量进行 (1) 中的计数, 最后将两者加起来返回即可.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int match(int *P,int n,int threhold){
    if(n==0)return 0;
    int cnt=0;
    if (n<=threhold){
        for(int i=n-1;i>=0;i--){
            if(P[i]==')')cnt++;
            if(P[i]=='(')cnt--;
        }
        return cnt;
    }
    par_do:
    	int cnt1=match(P,n/2);
    	int cnt2=match(P+n/2,n-n/2);
    return cnt1+cnt2;
}

仍然满足 , .

得到 .

3 Compute the Integral

采用分区间的方式计算积分 (1) 如果直接将

1
2
3
4
5
for (size_t i = 0; i < n; i++){
    double x = low + (high - low) * i / n; // 用的是每个小区间的左端点
    double dx = (high - low) / n;
    ans += f(x) * dx;
}

改成并行的 parallel_for:

1
2
3
4
5
6
auto cal = [&](size_t i) {
    double x = low + (high - low) * i / n;
    double dx = (high - low) / n;
    ans += f(x) * dx;
};
parallel_for(0, n, cal);

当 时, 两者的结果分别为

1
2
3
4
5
6
7
8
9
10
11
12
Total sum: 0.3740279091741
Warmup round running time: 1.875804
Round 1 running time: 2.020929
Round 2 running time: 2.056737
Round 3 running time: 1.742642
Average running time: 1.940102666667
Total sum: 0.1339990181118
Warmup round running time: 1.995736
Round 1 running time: 2.017125
Round 2 running time: 2.14402
Round 3 running time: 2.092346
Average running time: 2.084497

时间上并没有明显改变, 并且结果是错的.

原因很明显, 出现了 race, 如果两个进程同时执行 ans+=f(x)*dx, 这两个进程同时读到相同的 ans 并将自己的增量加到 ans 上, 由于没有互斥, 最终一个进程的结果会覆盖掉另一个进程的结果, 导致最终积分的结果仅仅是部分项的和. 还有就是这里每一项的计算 overhead 仍然不大, 因此 fork 的 overhead 还是处于支配地位.

  1. 解决上面问题最直接的方法是加锁.
1
2
3
4
5
6
7
8
#include<mutex>
std::mutex ans_mutex;
auto cal = [&](size_t i) {
    double x = low + (high - low) * i / n;
    double dx = (high - low) / n;
    std::lock_guard<std::mutex> lock(ans_mutex);
    ans += f(x) * dx;
};

这会破坏 parallelism, 下面是新的结果

1
2
3
4
5
6
7
8
9
10
11
12
Total sum: 0.3740279091741
Warmup round running time: 1.725927
Round 1 running time: 1.718653
Round 2 running time: 1.720002
Round 3 running time: 1.717201
Average running time: 1.718618666667
Total sum: 0.3740279091742
Warmup round running time: 10.593743
Round 1 running time: 10.475321
Round 2 running time: 10.407566
Round 3 running time: 10.49787
Average running time: 10.46025233333

时间变成了原来(串行)的五倍.

  1. 设计新的并行算法.

这就是计算 项和, 所以可以直接使用 reduce 算法, 同时增加粒度的控制:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
template <class Func>
double integral_reduce(const Func &f, size_t n, double low, double high, size_t threhold){
    auto rect = [&](size_t i) -> double{
        double x = low + (high - low) * i / n;
        double dx = (high - low) / n;
        return f(x) * dx;
    };
    if (n < threhold){
        double ans = 0;
        for (size_t i = 0; i < n; i++){
            ans += rect(i);
        }
        return ans;
    }
    else{
        double v1, v2;
        auto f1 = [&](){
            v1 = integral_reduce(f, n / 2, low, low + (high - low) * (n / 2) / n, threhold);
        };
        auto f2 = [&](){
            v2 = integral_reduce(f, n - n / 2, low + (high - low) * (n / 2) / n, high, threhold);
        };
        par_do(f1, f2);
        return v1 + v2;
    }
}

, 粒度为 , 8个线程时偶尔可以达到3秒内:

1
2
3
4
5
6
Total sum: 0.3740279120104
Warmup round running time: 2.975084
Round 1 running time: 2.995158
Round 2 running time: 2.970043
Round 3 running time: 2.956155
Average running time: 2.973785333333

粒度调细一些可能会更好.

  1. 再提高速度和精度!

上面的结果精度不够, 提升精度的一个很显然的方法就是用梯形代替矩形:

1
2
3
4
5
6
7
auto trape = [&](size_t i) -> double
{
    double l = low + (high - low) * i / n; 
    double r=low + (high - low) * (i+1) / n;
    double dx = (high - low) / n;
    return (f(l)+f(r)) * dx/2;
};

或者更简单的写法:

1
2
3
4
5
6
auto trape = [&](size_t i) -> double
{
    double l = low + (high - low) * i / n; 
    double dx = (high - low) / n;
    return (f(l)+f(l+dx)) * dx/2;
};

这样很容易就达到了要求的精度, 但是计算梯形的面积增加了 overhead, 时间变成原来的两倍了:

1
2
3
4
5
6
Total sum: 0.3740279123255
Warmup round running time: 7.468748
Round 1 running time: 7.875462
Round 2 running time: 8.723461
Round 3 running time: 7.854643
Average running time: 8.151188666667

很容易发现开销主要在除 之外的点都被算了两遍, 而且 dx 只需要计算一次. 同时由于开不了大小为 的数组存下所有的点(这样写还需要考虑区间的起始下标), 因此可以算完第 个梯形的面积后把下底传出作为第 个梯形的上底, 保证区间上的每个点只被计算一次, 如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
double dx = (high - low) / n;
auto trape = [&](double &ls, double r) -> double {
    double rs = f(r);
    double res = (ls + rs) * dx / 2;
    ls = rs;
    return res;
};
if (n < threhold) {
    double ans = 0, ls = f(low), r = low;
    for (size_t i = 0; i < n; i++) {
        r += dx;
        ans += trape(ls, r);
    }
    return ans;
}

最终 , 8 个线程, threhold=200 的结果如下:

1
2
3
4
5
6
7
8
9
10
11
12
Total sum: 0.3740279335276
Warmup round running time: 10.416636
Round 1 running time: 10.530778
Round 2 running time: 10.372824
Round 3 running time: 10.660854
Average running time: 10.52148533333
Total sum: 0.3740279123255
Warmup round running time: 2.727995
Round 1 running time: 2.793685
Round 2 running time: 2.812348
Round 3 running time: 2.778975
Average running time: 2.795002666667

第一个结果是串行计算, 第二个结果是并行计算, 两者都是按照梯形计算的.

为什么串行计算的结果偏大 ?

怎么搞到 2s 内 ?

REF :

[1] CS214 of UCR

[2] dropbox

Experiment 1: Parallel Eratosthenes Sieve (MPI)

串行的埃筛如下 (筛区间)

  1. (第一个未被标记的数)
  2. 当 时标记 中所有 的倍数转 3; 反之结束, 所有未被标记的数为素数.
  3. 下一个未被标记的数, 这一定是当前 的下一个素数.

复杂度为 .

为什么只需要标记 中所有 的倍数? 这意味着 中的合数已经被小于 的素数标记过了, 这是显然的, 因为对每一个 且 为合数, 一定有小于 的素因子, 否则 .

复杂度的证明:

参考 当数论遇上分析（3）——数论函数的加权平均、切比雪夫定理以及埃氏筛 - 知乎 (zhihu.com)

Mertens 定理 - 香蕉空间 (bananaspace.org)

黎曼－斯蒂尔杰斯积分有什么存在的意义吗？ - 知乎 (zhihu.com)

数论函数(1) - 知乎 (zhihu.com)

筛法 - OI Wiki (oi-wiki.org)

需要分别证明 Mertens 第一定理 和 Mertens 第二定理 .

这是典型的 RS 积分的应用, 同时 RS 积分和 Abel 变换又是一脉相承的. 

1

ok 先放一下, 先看代码.

MPI Something

1
2
3
4
5
6
7
int MPI_Bcast(
    void *buffer,          /* Addr of 1st element */
    int count,             /* # elements to broadcast */
    MPI_Datatype datatype, /* Type of elements */
    int root,              /* ID of root process */
    MPI_Comm comm)         /* Communicator */
MPI_Bcast(&k, 1, MPI_INT, 0, MPI_COMM_WORLD);

并行思路

很容易能发现标记的过程可以并行, 可以分块, 每个处理器分别负责不同区间的标记, 设有 个处理器:

  • 希望每个区间的大小为 或 保证负载均衡

  • 当 时, 直接均分即可

  • 当 时, 记 , 则 , 有 个长度为 的区间, 有 个长度为 的区间. 然后需要求出区间 的起始元素的下标, 这有很多组织方式, 需要找出较为表达式简单的那一种.

  • 区间 的第一个元素为 , 最后一个元素为 . 显然总的元素的个数为 同时 , 其中 , 注意到 那么 两不等式相加有 显然 即 上的整数仅有 两个 (也可能是一个), 那么就一定有 , 这就意味着每个区间的元素个数要么为 , 要么为 , 这就是我们所希望的.

  • 对于元素 而言, 它所属的区间为 .

设 个块为 , 块 的长度为且元素为 . 现在探究一下具体并行标记的过程.

从块首元素开始枚举, 找到第一个 的倍数标记, 接下来下标每次加 标记即可.

样例代码, 先不管.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <mpi.h>
#include <math.h>
#include <stdio.h>
#include "MyMPI.h"
#define MIN(a, b) ((a) < (b) ? (a) : (b))
int main(int argc, char *argv[])
{
    MPI_Init(&argc, &argv);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();
    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    if (argc != 2)
    {
        if (!id)
            printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }
    n = atoi(argv[1]);
    low_value = 2 + BLOCK_LOW(id, p, n - 1);
    high_value = 2 + BLOCK_HIGH(id, p, n - 1);
    size = BLOCK_SIZE(id, p, n - 1);
    proc0_size = (n - 1) / p;
    if ((2 + proc0_size) < (int)sqrt((double)n))
    {
        if (!id)
            printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }
    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 = 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])
                ; /*find the smallest unmarked*/
            prime = index + 2;
        }
        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);
    elapsed_time += MPI_Wtime();
    if (!id)
    {
        printf("%d primes are less than or equal to %d\n",
               global_count, n);
        printf("Total elapsed time: %10.6f\n", elapsed_time);
    }
    MPI_Finalize();
    return 0;
}

优化

  • 去掉偶数

  • 优化通信

  • 考虑 cache miss 重新组织循环

更多优化思路

  • bitset
  • 向量化?
  • 硬件上有无其他优化?
  • 想试试 wheel factorization, 但好像算是改算法了, 如果想在算法层面上做优化那么肯定就只能避免重复筛了, 反正肯定是比不过 .
  • 可能需要具体再了解了解 MPI

Wheel Factorization 也基于分块的思想.

我们先选取一个合数 并先筛区间 这时，对于区间 前 个区间中，每个区间的所有整数都是一个模 完全剩余系，对于 有 显然 为素数的一个必要条件是 因此我们在筛后面的块的时候只需要对满足 的数进行判断即可，这样的 有 个，此时我们已经可以直接写出复杂度了，复杂度为 现在要做的是设计合理的 使得 尽可能小，设 的质因数分解为 则 显然 越大 越小，且 与 无关，于是令 ，这样做预处理出 上的素数的时间最短。在 范围内，可以选择 这样已经能满足我们的要求了。此时

Parallel Floyd Algorithm

Floyd 算法计算 的每一对顶点的最短路径的长度, 算法需要三重循环. 设 表示 之间的最短路径长度, 则枚举所有的 , 更新 , 初始化 , 其中 为边 的边权.

MPI 点对点通信

2024-01-25 该篇文章被 Baoduo Xu 归为分类: Parallel Computing

以上
©
Too Young Too Naive 本站由Hexo驱动｜使用Hexo-theme-A4主题