學校課程/平行程式

數據分割 區間分塊小技巧

前言

運行平行程式之前，要將資料切成很多塊，每一塊之間獨立。從邏輯上看來有 row-based, column-based, matrix-based, tree-based 等平行方法。在許多方法中，最容易入手的就是 row-based。

假設有 $n$ 個元素，編號介於 $[0, n-1]$ 之間，要分成 $m$ 組。每一組得編號都要連續，以增加資料局部性 (data locality) 並且每一塊大小盡可能一樣，請問要怎麼分割才好？這樣的數學問題要做出來並不是件難事，但強迫症的人來說，在程式中零碎的判斷相當折騰。

數學分析

直觀

若有 $n$ 個元素分成 $m$ 組，每一個組至少有 $\lfloor n / m \rfloor$ 個元素，其中 $n \mod m$ 組會額外多一個元素。因此，對於第 $i$ 組滿足 $i < n \mod m$ 都會多一個元素。

1
2
3
4
5
6
7
8
9
10
void method0(int n, int m) {
    printf("In %s\n", __func__);
    int bsz = n/m, rm = n%m, sz;
    for (int i = 0, L = 0, R; i < m; i++) {
        sz = bsz + (i < rm);
        R = L + sz - 1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, (R - L + 1));
        L = R + 1;
    }
}

程序流

這方式屬於懶惰程序員，根據整數的性質會得到刻度，在計算第 $i$ 個刻度採用 $\lfloor \frac{in}{m} \rfloor$，自然而然地可以與下一個刻度形成一個區間。不幸地，若 $m > n$ 時，有些刻度會重疊，導致重複計算的情況發生。

1
2
3
4
5
6
7
void method1(int n, int m) {
    printf("In %s\n", __func__);
    for (int i = 0, L, R; i < m; i++) {
        L = i*n/m, R = (i+1)*n/m - 1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, (R - L + 1));
    }
}

ceil

從《具體數學》第三章節的介紹，得到幾個簡單的數學公式

$$\begin{align} n = \left \lceil \frac{n}{m} \right \rceil + \left \lceil \frac{n-1}{m} \right \rceil + \cdots + \left \lceil \frac{n-m-1}{m} \right \rceil \end{align}$$

第 $i$ 團 ($0 \le i < m$)，要處理 $\left \lceil \frac{n-i}{m} \right \rceil$ 個元素。

1
2
3
4
5
6
7
8
9
void method2(int n, int m) {
    printf("In %s\n", __func__);
    for (int i = 0, L = 0, R, sz; i < m; i++) {
        sz = (n + m-i-1) / m;
        R = L + sz-1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, sz);
        L = R + 1;
    }
}

floor

又或者使用 floor 表示法

$$\begin{align} n = \left \lfloor \frac{n}{m} \right \rfloor + \left \lfloor \frac{n+1}{m} \right \rfloor + \cdots + \left \lfloor \frac{n+m-1}{m} \right \rfloor \end{align}$$

第 $i$ 團 ($0 \le i < m$)，要處理 $\left \lceil \frac{n+i}{m} \right \rceil$ 個元素。

1
2
3
4
5
6
7
8
9
void method3(int n, int m) {
    printf("In %s\n", __func__);
    for (int i = 0, L = 0, R, sz; i < m; i++) {
        sz = (n + i) / m;
        R = L + sz-1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, sz);
        L = R + 1;
    }
}

補充流

不管怎樣，前 $i$ 就拿取 $\lceil \frac{n}{m} \rceil$，最後一組拿少一點。

1
2
3
4
5
6
7
8
9
void method4(int n, int m) {
    printf("In %s\n", __func__);
    int bsz = (n + m-1)/m, sz;
    for (int i = 0, L = 0, R; i < m; i++) {
        L = i*bsz, R = L + bsz - 1;
        if (R >= n)	R = n-1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, (R - L + 1));
    }
}

測試

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
#include <bits/stdc++.h>
using namespace std;
void method0(int n, int m) {
    printf("In %s\n", __func__);
    int bsz = n/m, rm = n%m, sz;
    for (int i = 0, L = 0, R; i < m; i++) {
        sz = bsz + (i < rm);
        R = L + sz - 1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, (R - L + 1));
        L = R + 1;
    }
}
void method1(int n, int m) {
    printf("In %s\n", __func__);
    for (int i = 0, L, R; i < m; i++) {
        L = i*n/m, R = (i+1)*n/m - 1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, (R - L + 1));
    }
}
void method2(int n, int m) {
    printf("In %s\n", __func__);
    for (int i = 0, L = 0, R, sz; i < m; i++) {
        sz = (n + m-i-1) / m;
        R = L + sz-1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, sz);
        L = R + 1;
    }
}
void method3(int n, int m) {
    printf("In %s\n", __func__);
    for (int i = 0, L = 0, R, sz; i < m; i++) {
        sz = (n + i) / m;
        R = L + sz-1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, sz);
        L = R + 1;
    }
}
void method4(int n, int m) {
    printf("In %s\n", __func__);
    int bsz = (n + m-1)/m, sz;
    for (int i = 0, L = 0, R; i < m; i++) {
        L = i*bsz, R = L + bsz - 1;
        if (R >= n)	R = n-1;
        printf("%d len([%d, %d]) = %d\n", i, L, R, (R - L + 1));
    }
}
int main() {
    void (*FUNC[])(int, int) = {method0, method1, method2, method3, method4};
    int n, m;
    while (scanf("%d %d", &n, &m) == 2) {
        for (int i = 0; i < sizeof(FUNC)/sizeof(FUNC[0]); i++)
            (*FUNC[i])(n, m);
    }
    return 0;
}
Read More +
學校課程/平行程式

批改娘 10085. Parallel Count (debug)

題目描述

這有一份由 pthread 撰寫的序列總和計算，假設不開任何的優化參數，在快取處理會有嚴重缺失。

main.c

輸入序列長度 $n$，計算 $m$ 次經由 $\text{key}, \; \text{key} + 1, \; \text{key} + 2, \; \cdots, \; \text{key} + m-1$ 加密的序列總和。這部份將不提供修改。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include <stdio.h>
#include <stdlib.h>
#include "utils.h"
int main() {
    int cases = 0;
    int n, m, key;
    scanf("%d %d %d", &n, &m, &key);
    for (int it = 0; it < m; it++) {
        int ret = run(n, key);
        printf("Case #%d: %d\n", ++cases, ret);
        key++;
    }
    return 0;
}

utils.h

計算工作交由 void f(int n, int key, int *p1, int *p2, int *p3, int *p4); 完成最後四個參數，將會由 4 個 thread 計算分別儲存在位址 p1, p2, p3, p4 中。

1
2
3
4
5
6
7
#ifndef __UTILS_H
#define __UTILS_H
void f(int n, int key, int *p1, int *p2, int *p3, int *p4);
int run(int n, int key);
#endif

sum.c

平行計算序列總和，特別注意到 void* subtask(void* argu) 中存取的記憶體位置。

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
#include <stdlib.h>
#include <pthread.h>
#include <stdint.h>
#include "utils.h"
typedef struct Argu {
    int *result;
    int L, R, key;
} Argu;
static inline uint32_t rotate_left(uint32_t x, uint32_t n) {
        return  (x << n) | (x >> (32-n));
}
static inline uint32_t encrypt(uint32_t m, uint32_t key) {
        return (rotate_left(m, key&31) + key)^key;
}
static void* subtask(void* argu) {
    Argu *ptr = (Argu *) argu;
    *(ptr->result) = 0;
    for (int i = ptr->L; i <= ptr->R; i++)
        *(ptr->result) += encrypt(i, ptr->key);
}
void f(int n, int key, int *p1, int *p2, int *p3, int *p4) {
    pthread_t threads[4];
    int *output[4] = {p1, p2, p3, p4};
    for (int i = 0; i < 4; i++) {
            Argu *argu = (Argu *) malloc(sizeof(Argu));
            argu->result = output[i];
            argu->L = i * (n / 4) + 1;
            argu->R = argu->L + n / 4 - 1;
            argu->key = key;
            if (i == 3) argu->R = n;
            pthread_create(&threads[i], NULL, subtask, argu);
    }
    for (int i = 0; i < 4; i++)
        pthread_join(threads[i], NULL);
}

job.c

你的工作要改善下方代碼的效能。

1
2
3
4
5
6
7
8
9
10
#include "utils.h"
int ret[128];
int run(int n, int key) {
    int sum = 0;
    f(n, key, ret, ret+1, ret+2, ret+3);
    for (int i = 0; i < 4; i++)
        sum += ret[i];
    return sum;
}

輸入格式

輸入只有一行三個整數 $n, \; m, \; \text{key}$。

輸出格式

輸出 $m$ 序列總和結果 (無視 overflow)。

範例輸入

1
10000000 10 514

範例輸出

1
2
3
4
5
6
7
8
9
10
Case #1: 1397862656
Case #2: 1970821632
Case #3: -1178356736
Case #4: 1113221120
Case #5: 1401409536
Case #6: 1977786368
Case #7: -1164427264
Case #8: 1145914243
Case #9: 645957382
Case #10: 1308383748

編譯環境

Makefile

1
2
3
4
5
6
CFLAG=-std=c99 -pthread
all: main.c sum.c job.c
    gcc $(CFLAG) main.c -c
    gcc $(CFLAG) sum.c -c
    gcc $(CFLAG) main.o sum.o job.c -o job

reference

Solution

由於茵可大神給予我參考的簡報 Mind Cache 中提到不少的快取實驗，其中一部份就是在 thread 使用上，於是就拿來出題測試測試一番。

通常講到數據局部性 (Data Locality) 都希望資料盡量集中，但一不小心會犯下錯誤，就是在多個 thread 儲存答案時，雖然不會共用同一個位址，但相鄰的位址在另一個核 (core) 使用，由於彼此之間相互修改，兩個相鄰位址若放在同一個 cache line，dirty bit 勢必要讓他們從 cache line 掉到 L1 -> L2 -> L3，進行資料同步。

目前 cache line 普遍上設計大小為 64 Bytes，相當於間隔 16 個 4 bytes Integer 的差距，所以儲存答案間隔遠一點會比近一點好。當然，編譯器優化參數開到 O2 以上時，似乎就直接先存放到暫存器裡頭，不會不斷地存取 global memory 部份，如此一來，就不會發生 cache miss 問題。實驗結果效能可以差到四倍左右。

1
2
3
4
5
6
7
8
9
#include "utils.h"
 
int ret[128];
int run(int n, int key) {
    int sum = 0;
    f(n, key, ret+0, ret+16, ret+32, ret+48);
    sum = ret[0] + ret[16] + ret[32] + ret[48];
    return sum;
}
Read More +
學校課程/平行程式

批改娘 10084. Prefix Sum (pthread)

題目描述

給予序列 $A \left[ 1 \cdots n \right]$，請計算前綴和序列 $S \left[ 1 \cdots n \right]$，其中 $$S \left[ i \right] = \sum_{k=1}^{i} A \left[ k \right]$$

為了檢測方便，移除輸入和輸出時間，序列 $A$ 由一個簡單加密 $\textit{key}$ 得到序列 $$A \left[ i \right] = \textit{encrypt}(i, \textit{key})$$以及輸出部分直接呼叫$\textit{output}(\textit{S}, n)$。注意$S\left[ 0 \right] = 0$，儲存答案的範圍為$S\left[ 1 \cdots n \right]$。

utils.h

可以直接 #include "utils.h" 在你的程式中。

1
2
3
4
5
6
7
8
9
10
11
#ifndef _UTILS_H
#define _UTILS_H
#include <stdint.h>
static inline uint32_t rotate_left(uint32_t x, uint32_t n) {
    return  (x << n) | (x >> (32-n));
}
static inline uint32_t encrypt(uint32_t m, uint32_t key) {
    return (rotate_left(m, key&31) + key)^key;
}
void output(uint32_t presum[], int n);
#endif

prefixsum-seq.c

請修改這份程序。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include "utils.h"
#define MAXN 10000005
#define MAX_THREAD 4
uint32_t prefix_sum[MAXN];
int main() {
    int n;
    uint32_t key;
    while (scanf("%d %" PRIu32, &n, &key) == 2) {
        uint32_t sum = 0;
        for (int i = 1; i <= n; i++) {
            sum += encrypt(i, key);
            prefix_sum[i] = sum;
        }
        output(prefix_sum, n);
    }
    return 0;
}

secret.c (測試用)

實際情況會用助教寫好的平行方式進行計算且雜湊公式不同。

1
2
3
4
5
6
7
8
9
10
11
12
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include "utils.h"
void output(uint32_t presum[], int n) {
    uint32_t hash = 0;
    for (int i = 1; i <= n; i++)
        hash += presum[i] * i;
    printf("%" PRIu32 "\n", hash);
}

輸入格式

輸入有多組測資，每組測資一行兩個整數 $n$, $\textit{key}$，分別為序列長度 $n$，加密金鑰 $\textit{key}$。

  • $1 \le n \le 10^7$
  • $0 \le key < 2^{32}$

輸出格式

對於每一組測資呼叫 output(uint32_t presum[], int n)，隨後會輸出一個數值。

範例輸入

1
2
3 2
10 5

範例輸出

1
2
100
54560

範例解釋

  • $(n, \textit{key})=(3, 2)$，得 $A \left[ 1 \cdots 3\right] = \left[ 4, 8, 12 \right]$，$S \left[ 1 \cdots 3\right] = \left[ 4, 12, 24 \right]$，$\text{hash} = 4 + 12 \times 2 + 24 \times 3 = 100$

編譯方式

1
gcc -std=c99 -O2 -pthread prefixsum-seq.c secret.c -o prefixsum-seq

測試主機資訊

推薦使用 4 個執行緒解決此問題，平行效能接近 2 倍改進。若使用過多的執行緒，將因為要排到另一個處理器上運行而變慢。

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
processor   : 0
vendor_id   : GenuineIntel
cpu family  : 6
model       : 62
model name  : Intel(R) Xeon(R) CPU E5-2620 v2 @ 2.10GHz
stepping    : 4
microcode   : 0x415
cpu MHz     : 1200.000
cache size  : 15360 KB
physical id : 0
siblings    : 12
core id     : 0
cpu cores   : 6
apicid      : 0
initial apicid  : 0
fpu     : yes
fpu_exception   : yes
cpuid level : 13
wp      : yes
processor   : 1
vendor_id   : GenuineIntel
cpu family  : 6
model       : 62
model name  : Intel(R) Xeon(R) CPU E5-2620 v2 @ 2.10GHz
stepping    : 4
microcode   : 0x415
cpu MHz     : 1200.000
cache size  : 15360 KB
physical id : 0
siblings    : 12
core id     : 1
cpu cores   : 6
apicid      : 2
initial apicid  : 2
fpu     : yes
fpu_exception   : yes
cpuid level : 13
wp      : yes

Solution

利用 $P$ 的 thread，分開計算區間 $[1, N/P], \; [N/P+1, 2N/P], \cdots$ 的前綴和，由於前綴和只有利用加法計算，加法具有結合律，隨後傳遞前 $i$ 段的和，可以循序 $\mathcal{O}(P)$ 完成，或者利用樹形平行 (tree-based) 的方式在 $\mathcal{O}(\log P)$，在實務上直接循序即可。

由於採用分治方式，需要平行兩次，時間複雜度為 $\mathcal{O}(2 \frac{N}{P} + P)$。

  • 當 $P = 3$ 時，只獲得 1.5 倍的加速
  • 當 $P = 4$ 時，獲得 2 倍加速
  • 當 $P > 4$ 時，要將記憶體的部分傳到另一顆 CPU 上，假設搬運時間是線性正比，計算複雜度也是線性時間，那麼搬運時間不如在同一顆 CPU 上運行。
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
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <pthread.h>
#include "utils.h"
 
typedef struct Argu {
    int l, r;
    uint32_t *psum;
    uint32_t key;
} Argu;
int min(int x, int y) {
    return x < y ? x : y;
}
void* subtask1(void *argu) {
    Argu data = *((Argu *) argu);
    int l = data.l, r = data.r;
    uint32_t *psum = data.psum; 
    uint32_t sum = 0, key = data.key;
    for (int i = l, j = 0; i <= r; i++, j++) {
        sum += encrypt(i, key);
        psum[j] = sum;
    }
    free(argu);
}
void* subtask2(void *argu) {
    Argu data = *((Argu *) argu);
    int l = data.l, r = data.r;
    uint32_t *psum = data.psum; 
    uint32_t base = data.key;
    for (int i = l, j = 0; i <= r; i++, j++)
        psum[j] += base;
    free(argu);
}
#define MAXN 10000005
#define MAX_THREAD 4
uint32_t prefix_sum[MAXN];
int main() {
    int n;
    uint32_t key;
    while (scanf("%d %" PRIu32, &n, &key) == 2) {
        int BLOCK = (n+MAX_THREAD-1) / MAX_THREAD, m = 0;
        pthread_t threads[MAX_THREAD];
        for (int i = 1; i <= n; ) {
            Argu *argu = (Argu *) malloc(sizeof(Argu));
            argu->l = i, argu->r = min(n, i + BLOCK - 1);
            argu->psum = prefix_sum + i, argu->key = key;
            i += BLOCK;
            pthread_create(&threads[m], NULL, subtask1, argu), m++;
        }
        for (int i = 0; i < m; i++)
            pthread_join(threads[i], NULL);
 
        m = 0;
        uint32_t block_sum = 0;
        for (int i = 1; i <= n; i) {
            uint32_t tmp = block_sum;
            block_sum += prefix_sum[min(n, i+BLOCK-1)];
            Argu *argu = (Argu *) malloc(sizeof(Argu));
            argu->l = i, argu->r = min(n, i + BLOCK - 1);
            argu->psum = prefix_sum + i, argu->key = tmp;
            i += BLOCK;
            pthread_create(&threads[m], NULL, subtask2, argu), m++;
        }
        for (int i = 0; i < m; i++)
            pthread_join(threads[i], NULL);
        output(prefix_sum, n);
    }
    return 0;
}
Read More +