2020q1 Homework2 (fibdrv)

contributed by < OscarShiang >

tags: `linux2020`

測試環境

$ cat /etc/os-release 
NAME="Ubuntu"
VERSION="18.04.4 LTS (Bionic Beaver)"
ID=ubuntu
ID_LIKE=debian
PRETTY_NAME="Ubuntu 18.04.4 LTS"
VERSION_ID="18.04"
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
VERSION_CODENAME=bionic
UBUNTU_CODENAME=bionic

$ cat /proc/version
Linux version 5.3.0-40-generic (buildd@lcy01-amd64-024) 
(gcc version 7.4.0 (Ubuntu 7.4.0-1ubuntu1~18.04.1)) 
#32~18.04.1-Ubuntu SMP Mon Feb 3 14:05:59 UTC 2020

作業要求

研讀上述 Linux 效能分析的提示描述，在自己的實體電腦運作 GNU/Linux，做好必要的設定和準備工作
$\to$ 從中也該理解為何不希望在虛擬機器中進行實驗;
研讀上述費氏數列相關材料 (包含論文)，摘錄關鍵手法，並思考 clz / ctz 一類的指令對 Fibonacci 數運算的幫助。請列出關鍵程式碼並解說
複習 C 語言數值系統和 bitwise operation，思考 Fibonacci 數快速計算演算法的實作中如何減少乘法運算的成本;
lsmod 的輸出結果有一欄名為 Used by，這是 "each module's use count and a list of referring modules"，但如何實作出來呢？模組間的相依性和實際使用次數 (reference counting) 在 Linux 核心如何追蹤呢？
注意到 fibdrv.c 存在著 DEFINE_MUTEX, mutex_trylock, mutex_init, mutex_unlock, mutex_destroy 等字樣，什麼場景中會需要呢？撰寫多執行緒的 userspace 程式來測試，觀察 Linux 核心模組若沒用到 mutex，到底會發生什麼問題
修正 Fibonacci 數計算的正確性，改善 fibdrv 計算 Fibinacci 數列的執行效率，過程中需要量化執行時間，可在 Linux 核心模組和使用層級去測量

改善 Fibonacci 的計算效率

Fibonacci 數的定義如下：

F_{0} = 1, F_{1} = 1 F_{n} = F_{n - 1} + F_{n - 2} (n \geq 2)

而在程式中的實作為

int F(int n) 
{
    if (n == 0 || n == 1)
        return 1;
    return F(n - 1) + F(n - 2);
}

假設計算 F(6) 的時候，最後程式呼叫的結果就會變成下圖:

在該展開圖裡，F(4), F(3), F(2) 在遞迴的過程中就會被重複計算，很顯然這樣的計算效率非常的差。根據 Vorobev 等式，我們可以將 Fibonacci 數的計算方式整理成下列形式:

[\begin{matrix} F_{n} \\ F_{n - 1} \end{matrix}] = [\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}] \times [\begin{matrix} F_{n - 1} \\ F_{n - 2} \end{matrix}]

而

F_{2 n + 1}

與

F_{2 n}

則變成

\begin{aligned} [\begin{array}{c} F_{2 n + 1} \\ F_{2 n} \end{array}] & = {[\begin{array}{c} 1 & 1 \\ 1 & 0 \end{array}]}^{2 n} \times [\begin{array}{c} F_{1} \\ F_{0} \end{array}] \\ = {[\begin{array}{c} 1 & 1 \\ 1 & 0 \end{array}]}^{n} \times {[\begin{array}{c} 1 & 1 \\ 1 & 0 \end{array}]}^{n} \times [\begin{array}{c} F_{1} \\ F_{0} \end{array}] \\ = [\begin{array}{c} F_{n + 1} & F_{n} \\ F n & F_{n - 1} \end{array}] \times [\begin{array}{c} F_{n + 1} & F_{n} \\ F n & F_{n - 1} \end{array}] \times [\begin{array}{c} 1 \\ 0 \end{array}] \\ = [\begin{array}{c} F_{n + 1}^{2} + F_{n}^{2} \\ F_{n} (2 F_{n + 1} - F_{n}) \end{array}] \end{aligned}

所以

F_{2 n + 1}

與

F_{2 n}

就可以表示為

F_{2 n + 1} = F_{n + 1}^{2} + F_{n}^{2} F_{2 n} = F_{n} (2 F_{n + 1} - F_{n})

只要我們有

F_{0} = 0, F_{1} = 1, F_{2} = 1

，這樣的整理可以避免逐項計算 Fibonacci 數所花費的時間:

同時因為電腦數值系統的特性，可以將其轉換為利用數值的 bit 來進行判斷：

遇到 0 : 計算
$F (2 n + 1)$ 與
$F (2 n)$
遇到 1 : 計算
$F (2 n + 1)$ 與
$F (2 n)$ 後再計算
$F (2 n + 2)$

使用 fibdrv 進行實驗

為了在 fibdrv 這個核心模組裡面測量使用原本的遞迴方式計算 fibonacci 與使用 Vorobev 等式計算 fibonacci 的時間差異，並進行分析

以下使用下列 fast doubling 的函式進行測試

static long long fib_doubling(long long k)
{
    if (k == 0)
        return 0;

    long long t0 = 1, t1 = 1, t3 = 1, t4 = 0;
    long long i = 1;
    while (i < k) {
        if ((i << 1) <= k) {
            t4 = t1 * t1 + t0 * t0; 
            t3 = t0 * (2 * t1 - t0);
            t0 = t3; 
            t1 = t4; 
            i <<= 1;
        } else {
            t0 = t3; 
            t3 = t4; 
            t4 = t0 + t4; 
            i++;
        }
    }   
    return t3; 
}

根據 The Linux Kernel Documentation 中所描述 ktime 的 API
我將 fib_read 這個函式改寫

/* calculate the fibonacci number at given offset */
static ssize_t fib_read(struct file *file,
                        char *buf,
                        size_t size,
                        loff_t *offset)
{
    // record the execution time by calling ktime_get()
    ktime_t kt = ktime_get();
    ssize_t result = fib_sequence(*offset);
    duration = ktime_sub(ktime_get(), kt);
    
    ktime_t kt2 = ktime_get();
    fib_doubling(*offset);
    kt2 = ktime_sub(ktime_get(), kt2);
    
    printk(KERN_INFO "%lld %lld", ktime_to_ns(kt), ktime_to_ns(kt2); 
    return result;
}

在執行 User Space 的 client 程式後利用終端機指令 dmesg 查看 fibdrv 的輸出

...
[ 1192.078203] 276 98
[ 1192.078210] 243 100
[ 1192.078217] 259 104
[ 1192.078223] 244 129
[ 1192.078230] 274 139
[ 1192.078237] 230 123
[ 1192.078244] 297 116
[ 1192.078251] 296 144
[ 1192.078258] 295 151
[ 1192.078265] 282 147
[ 1192.078272] 290 138
[ 1192.078279] 273 144
[ 1192.078286] 309 124
[ 1192.078293] 341 151
[ 1192.078300] 285 149
[ 1192.078307] 272 151

將這些實驗的數據以 gnuplot 繪製成圖表

在圖表中可以看到，兩者之間的差異在第 30 項前的差異並不大，但在第 30 項後差異持續擴大

根據 task switch 的描述：

In computing, a context switch is the process of storing the state of a process or thread, so that it can be restored and resume execution at a later point.

所以在進行上述實驗的時候就會導致 tail recursive 的花費時間分佈有抖動的情況發生

為了避免這樣的情況，我將函式改到 User Space ，並以 time.h 中的 struct timespec 計時

    for (int n = 0; n <= 100; n++) {
        long long fib;
        struct timespec ts, ts1;

        // test with adding version
        clock_gettime(CLOCK_REALTIME, &ts);
        fib = fib_adding(n);
        clock_gettime(CLOCK_REALTIME, &ts1);

        printf("time cost: %ld", ts1.tv_nsec - ts.tv_nsec);

        // test with doubling version
        clock_gettime(CLOCK_REALTIME, &ts);
        fib = fib_doubling(n);
        clock_gettime(CLOCK_REALTIME, &ts1);

        printf(" %ld\n", ts1.tv_nsec - ts.tv_nsec);
    }

並透過 taskset 將行程綁定在第 0 個 CPU 上執行

$ taskset 0x1 ./fib

最後將實驗結果繪製成下圖

從這張圖中可以看到，時間花費數值跳動的情況已經改善許多

透過設定 /boot/grub/grub.cfg 更改 bootloader 的參數，將開機指令 linux 後面加上 isolcpus=0 將第 0 個 CPU 排除在開機程序之外

重新開機後，確認第 0 個 CPU 已經被排除成功

利用 taskset 設定測試程式執行於第 0 個 CPU 上，測試結果如下

可以發現在限定使用特定 CPU 來進行實驗時，實驗所造成的誤差與干擾與在完全不調整的情況下相差非常多

為了避免單一實驗之極端值所造成的實驗誤差，我以 Python 腳本將實驗的過程與資料的整理自動化

#!/usr/bin/env python3

import os 
import pandas as pd
import math

PATH = '/home/oscar/Desktop/statistic'

times = 10000

adding = pd.DataFrame(columns = range(times))
doubling = pd.DataFrame(columns = range(times))

# executing experiment
for i in range(times):
	os.system('taskset 0x1 ./fib > ' + PATH + '/tmp/data' + str(i))

for i in range(times):
	f = open(PATH + '/tmp/data' + str(i), 'r')
	content = f.read()
	raw = content.split('\n')
	raw.remove('')

	add = []
	double = []
	for j in range(len(raw)):
		numbers = raw[j].split(' ')
		add.append(int(numbers[0]))
		double.append(int(numbers[1]))
	adding[i] = add;
	doubling[i] = double
	print('data' + str(i) + ' completed')

print(adding)
print('-' * 30)
print(doubling)

# processing the data
mean_add, std_add = adding.mean(axis = 1), adding.std(axis = 1)
mean_double, std_double = doubling.mean(axis = 1), doubling.std(axis = 1)

add = []
double = []

print('processing adding data...')
for i in range(len(adding.index)):
	sum = 0
	cnt = 0
	for j in adding.iloc[i]:
		if abs(j - mean_add[i]) < 2 * std_add[i]:
			sum = sum + j;
			cnt = cnt + 1;
	sum = sum / cnt;
	add.append(sum);

print('processing doubling data...')
for i in range(len(doubling.index)):
	sum = 0
	cnt = 0
	for j in doubling.iloc[i]:
		if abs(j - mean_double[i]) < 2 * std_double[i]:
			sum = sum + j;
			cnt = cnt + 1;
	sum = sum / cnt;
	double.append(sum);

out = []
for i, j, k in zip(range(len(add)), add, double):
	out.append('{} {} {}'.format(i, j, k))

# output the result
f = open(PATH + '/runtime_result', 'w')
print('file write: {}'.format(f.write('\n'.join(out))))
f.close()

重複實驗

10^{4}

次後利用 95% 的信賴區間統計資料，得到下圖

使用硬體手段加速

但不是所有的整數都會將其記憶體的 32 個位子填滿，所有我們可以利用 clz 計算 leading zeros 的個數，讓我們可以直接從數字的有效位元開始計算，而不需要一直判斷不需要的 leading zeros

根據 GCC 內建函式的描述，在 GCC 中就有直接可以使用的 clz / ctz 功能，所以我們可以在程式中透過 macro 將 clz 引進程式中

#define clz(s) __builtin_clz(s)

並利用 clz 先將 leading zeros 計算出來，作為 for 迴圈的中止點

根據 GCC 內建函式對於 __builtin_clz 回傳值的說明

Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.

所以為避免發生 undefined value 產生，在 n = 0 時直接回傳 0

因為 n = 0 的情況在計算的時候只會發生 2 次，所以利用 unlikely 巨集 (在 Linux 核心原始程式碼也存在) 提示 compiler 產生對 CPU 分支預測友善的程式碼。

將 unlikely 加入程式中

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

並將 n = 0 的判斷加入 unlikely 的標示

long long fib_doubling_with_clz(int n)
{
    if (unlikely(!n))
        return 0;

    long long a = 0, b = 1;
    for (int i = 32 - clz(n); i >= 0; i--) {
        long t1 = a * (2 * b - a); 
        long t2 = b * b + a * a;
        a = t1; 
        b = t2; 
        if (n & (1 << i)) {
            t1 = a + b;
            a = b;
            b = t1; 
        }
    }   
    return a;
}

與前面的實驗相同，本程式一樣使用 taskset 綁定來進行測量，而實驗結果如下圖所示

在上圖的實驗中，有 41 項 Fibonacci 數的計算是使用 clz 改善後的版本較佳，所以根據實驗結果可以知道 clz 的確能加速 fast doubling 的計算

透過重複實驗降低實驗誤差與極端值後結果呈現下圖

計算第 92 項之後的 Fibonacci 數

在進行 make check 時，會發現測試在第 93 項時會發生錯誤

利用測試程式輸出計算結果後發現在第 92 項之後出現 Overflow 的問題

...
90 2880067194370816120
91 4660046610375530309
92 7540113804746346429
93 -6246583658587674878
94 1293530146158671551
95 -4953053512429003327
96 -3659523366270331776
97 -8612576878699335103
98 6174643828739884737
99 -2437933049959450366
100 3736710778780434371

因為上述的情況發生在使用 long long 時，表示使用內建的整數型態以無法符合我們的要求，所以需要改以自定義的 struct 來改進結構

在 struct 的設計與運算實作上我參考〈Bignums Alrithmetic〉進行設計，使用 usigned int * 來存取一個動態長度的陣列

typedef struct _bign {
    unsigned int *d;
    unsigned int size;
} bn;

對照 quiz4，你可看到透過 small buffer optimization 手法實作的 vector，也能作為上述動態長度陣列用途，這也是隨堂測驗的規劃。
:notes: jserv

接著依序將加法、減法與陣列的調整操作實作出來

加法的部分，我參考〈Bignums Alrithmetic〉所提及的做法，利用 uint64_t 保留每個元素相加之後的進位。最後利用 15 ~ 20 行的 for 迴圈確認陣列的實際長度並將陣列的長度進行調整




























void bn_add(bn *a, bn b) 
{
    uint32_t len = max(a->size, b.size) + 1;
    uint32_t num[len];
    memset(num, 0, sizeof(num));

    uint64_t l = 0;
    for (int i = 0; i < len; i++) {
        uint32_t A = (i < a->size) ? a->d[i] : 0;
        uint32_t B = (i < b.size) ? b.d[i] : 0;
        l += (uint64_t)A + B;
        num[i] = l;
        l >>= 32; 
    }   
    for (int i = len - 1; i >= 0; i--) {
        if (!num[i])
            len--;
        else
            break;
    }   

    // adjustment array size
    if (len > a->size) {
        a->d = realloc(a->d, len * sizeof(int));
        a->size = len;
    }   
    memcpy(a->d, num, sizeof(int) * len);
}

但因為在 GCC 裡沒有對應 64 bit 以上 int 型態位元轉換的實作，所以我參考 How to convert a 128-bit integer to a decimal ascii string in C? 討論中，將大數轉換成十進位的字串形式，方便之後輸出運算的結果



































char *bn_todec(bn a) {
    char str[8 * sizeof(int) * a.size / 3 + 2];
    unsigned int n[a.size + 1];

    memset(str, '0', sizeof(str) - 1);
    str[sizeof(str) - 1] = '\0';

    memcpy(n, a.d, sizeof(int) * a.size);
    for (int i = 0; i < 8 * sizeof(int) * a.size; i ++) {
        int carry;

        carry = (n[a.size - 1] >= 0x80000000);

        for (int j = a.size - 1; j >= 0; j--)
            n[j] = ((n[j] << 1) & 0xffffffff) + ((j - 1) >= 0 ? (n[j - 1] >= 0x80000000) : 0);

        for (int j = sizeof(str) - 2; j >= 0; j--) {
            str[j] += str[j] - '0' + carry;
            carry = (str[j] > '9');
            if (carry)
                str[j] -= 10;
        }
    }

    // search for numbers
    int i = 0;
    while (i < sizeof(str) - 2 && str[i] == '0')
        i++;

    // passing string back
    char *p = malloc(sizeof(str) - i);
    memcpy(p, str + i, sizeof(str) - i);

    return p;
}

利用 bitwise 的運算將十進位的每位數字算出來後存進 str ，最後利用 27 行的迴圈計算 str 中未使用的數量，最後利用 malloc 取得一塊記憶體空間，並將 str 中有效的元素利用 memcpy 將十進位的字串指標 p 回傳回去

因為使用 fast doubling 計算 fibonacci 數的做法有使用到乘法，所以我使用 Karatsuba algorithm 來實作大數的乘法

因為 Karatsuba 用到 bitwise 的左、右位移，所以依需求實作出大數位元的位移

void bn_rshift(bn *out, int wid)
{
    int data = wid >> 5;  // devide by 32
    wid &= 31;            // mod 31

    // move int element across array
    if (data) {
        int tmp[data];
        memcpy(tmp, &out->d[data], sizeof(int) * (out->size - data));
        memset(out->d, 0, sizeof(int) * out->size);
        memcpy(out->d, tmp, sizeof(int) * (out->size - data));
        out->size -= data;
        out->d = realloc(out->d, sizeof(int) * out->size);
    }
    if (wid) {
        // move bit in the same array
        if (out->size == 1)
            out->d[0] >>= wid;
        // mvoe bit across array
        else {
            out->d[0] >>= wid;
            int mask = set_mask(wid);
            for (int i = 1; i < out->size; i++) {
                int part = out->d[i] & mask;
                if (part)
                    out->d[i - 1] |= (part << (32 - wid));
                out->d[i] >>= wid;
            }
        }
    }
    // clean up the redundant array
    for (int i = out->size - 1; i >= 1; i--) {
        if (!out->d[i])
            out->size--;
        else
            break;
    }
    out->d = realloc(out->d, out->size * sizeof(int));
}

右位移在實作上先處理超過 32 位元的位移，利用 memcpy 將整個陣列進行位移，當處理完超過 32 位元以上的位移後，處理小於 32 位元的位移，利用 mask 擷取跨陣列元素位移的位元，並將其利用 |= 放在對應的位置上，完成整個大數的位移。

void bn_lshift(bn *out, int wid)
{
    int data = wid / 32;
    wid %= 32;
    if (wid) {
        for (int i = out->size - 1; i >= 0; i--) {
            int zeros = clz(out->d[i]);
            if (wid <= zeros) {
                out->d[i] <<= wid;
            } else {
                // for fast lhs element -> realloc new
                if (i == out->size - 1) {
                    out->d = realloc(out->d, (out->size + 1) * sizeof(int));
                    memset(out->d + out->size, 0, sizeof(int));
                    ++out->size;
                }
                out->d[i] <<= zeros;
                int mask = set_mask(wid - zeros) << (32 - wid + zeros);
                int tmp = (mask & out->d[i]) >> (32 - wid + zeros);
                out->d[i + 1] |= tmp;
                out->d[i] <<= (wid - zeros);
            }
        }
    }
    if (data) {
        uint32_t tmp[out->size];
        memcpy(tmp, out->d, out->size * sizeof(int));
        out->d = realloc(out->d, (out->size + data) * sizeof(int));
        memset(out->d, 0, (out->size + data) * sizeof(int));
        memcpy(out->d + data, tmp, out->size * sizeof(int));
        out->size += data;
    }
}

左位移在實作上則與右位移有些不同：左位移首先處理小於 32 位元的位移，並利用 clz 檢查是否會導致 msb 溢位的問題產生，若會產生問題，則調整陣列大小，並將 msb 移動到新的陣列元素上
接著再進行超過 32 位元的位移，一樣利用 memcpy 進行陣列元素長度的移動，最後將陣列大小寫入回 out->size 完成位移

完成左右位移的建構之後就可以實作出 Karatsuba 的乘法了

最後將 bn 的結構引進測試程式中，並改寫 fibonacci 的計算函式，測試結果如下：

...
n = 87, 679891637638612258
n = 88, 1100087778366101931
n = 89, 1779979416004714189
n = 90, 2880067194370816120
n = 91, 4660046610375530309
n = 92, 7540113804746346429
n = 93, 12200160415121876738
n = 94, 19740274219868223167
n = 95, 31940434634990099905
n = 96, 51680708854858323072
n = 97, 83621143489848422977
n = 98, 135301852344706746049
n = 99, 218922995834555169026
n = 100, 354224848179261915075

與 The first 300 Fibonacci numbers 的資料比對後確認計算結果無誤，並將其引進 fibdrv 中

將 fib_sequence 進行改寫

static char *fib_sequence(int k)
{
    bn out;
    bn_init(&out);

    fib_doubling(&out, k);

    char *str = bn_todec(out);
    bn_free(&out);
    return str;
}

因為最後的計算結果會以字串的方式回傳回來，所以將 fib_read 改寫為使用 linux/uaccess.h 所提供的 copy_to_user 函式

static ssize_t fib_read(struct file *file,
                        char *buf,
                        size_t size,
                        loff_t *offset)
{
    char *str = fib_sequence(*offset);
    copy_to_user(buf, str, strlen(str));
    return 0;
}

但因為在 The Linux Kernel API 中並沒有使用者層級常用的 malloc, free, memcpy 等函式，所以將 bign.c 加上 macro 改寫原先使用 stdlib.h 所進行的記憶體管理

可比照 memory.h，對記憶體操作函式包裝，你可以提供兩種版本，一個是對應到 <stdlib.h>，另一個則是 kmalloc, kzalloc 一類 Linux 核心函式。
:notes: jserv

好的，已新增 memory.h 來包裝記憶體管理函式
Oscar[color= orange]

在新增 memory.h 後，嘗試將 bign.c, bign.c 與 memory.h 在編譯過程中與 fibdrv.ko 連結起來，卻得到下列結果

...
ERROR: "bn_copy" [/home/oscar/linux2020/fibdrv/fibdrv.ko] undefined!
ERROR: "bn_free" [/home/oscar/linux2020/fibdrv/fibdrv.ko] undefined!
ERROR: "bn_assign" [/home/oscar/linux2020/fibdrv/fibdrv.ko] undefined!
ERROR: "bn_init" [/home/oscar/linux2020/fibdrv/fibdrv.ko] undefined!
WARNING: modpost: missing MODULE_LICENSE() in /home/oscar/linux2020/fibdrv/fibdrv.o

直接將 bign.c 中的所有函式定義直接加入 fibdrv.c 中進行編譯，但在執行 client 程式時，在計算 Fibonacci 數時電腦當機。

利用 Valgrind 觀察使用者層級的測試程式，推測是記憶體管理問題導致問題產生

2020q1 Homework2 (fibdrv)

tags: linux2020

測試環境

作業要求

改善 Fibonacci 的計算效率

使用 fibdrv 進行實驗

使用硬體手段加速

計算第 92 項之後的 Fibonacci 數

參考資料

Read more

Banana Pi R4 - Debricking Guide

Interesting Facts of Pthread

OS HW1: Benchmark Your Computer Black Box

程式設計(二) Lab 3 講解

tags: `linux2020`