平行程式設計(1):HW1: PI 蒙地卡羅演算法

HW1要使用蒙地卡羅演算法計算PI，其過程非常簡單，就像射鏢靶一樣隨機打點在座標上，而原面積相對於整個面積的比例可以看成是射中靶的機率，如果投擲的鏢非常多，那麼落在圓中的數量比上全部的數量應該要等於圓面積比上整個面積，如此一來可以計算出圓面積就可以計算出PI了。

而程式碼非常簡單，我們用rand()生成隨機座標然後把結果加總起來就好，但是因為這是平行程式所以我們要改用thread來加速，而rand()在thread中會出現問題，這等等我們會看到。

所以寫了一個普通版本的PI.exe，1百萬次的迭代執行時間大概落在 0.3S，然後馬上寫了thread版本結果執行時間高達3S差了10倍！開thread竟然還變慢！？會不會是工作站的電腦是單核心

呢？輸入cat /proc/cpuinfo之後發現工作站的電腦高達4核心！所以是程式的問題。

問題發現其實是在critical section放在每一個迭代中，所以浪費非常多的時間讓thread在wait lock，所以馬上把加總"toss_in_circle"這個變數改成local版本，然後我們在thread function結束時做一次加總到global_toss_in_circle的動作，因此critical section改到這裡就可以了。

本來以為這樣就沒問題了，但是花費時間從3變成0.5S，雖然變快，但是還是比沒用thread的版本慢，慢慢測試發現時間會一直跟著輸入toss數目增加，1千萬時PI_thread.exe 要花5S ，1億次時PI_thread.exe 要花50S ，所以這個花費時間果然還是在loop中，那是哪個地方會讓開啟thread時反而更費時呢？

沒錯就是rand()!查了一下rand()，這是一個non thread safe的posix function，原因是他用了static的變數，看一下他的definition在stdlib.h

static long holdrand = 1L;
...
int rand() {
  return (((holdrand = holdrand * 214013L + 2531011L) >> 16) & 0x7fff);
}

holdrand就是那個static的變數，他在使用srand(seed)時被定義

void srand(unsigned int seed) {
   holdrand = (long) seed;
}

所以這就是為什麼在thread中會更費時的原因，原因就是每個thread呼叫rand()的時候holdrand是共用的，當有其中一個thread寫入時，另一個process的thread的cache就要重新搬入記憶體，這個事費時的主要原因。

所以當我只有使用一個thread的時候花費時間跟PI.exe 1000000一樣幾乎是0.3S，但是一旦超過一個thread就會增加時間到0.5，2個thread跟3個thread跟4個thread時間一樣是0.5S所以重點不是thread多更費時，而是這個static variable讓process沒辦法快速執行。

於是應該要使用rand_r()這個thread safe的function來取得亂數，它的定義就沒有static的變數了

int
    rand_r (unsigned int *seed)
    {
      unsigned int next = *seed;
      int result;

      next *= 1103515245;
      next += 12345;
      result = (unsigned int) (next / 65536) % 2048;

      next *= 1103515245;
      next += 12345;
      result <<= 10;
      result ^= (unsigned int) (next / 65536) % 1024;

      next *= 1103515245;
      next += 12345;
      result <<= 10;
      result ^= (unsigned int) (next / 65536) % 1024;

      *seed = next;

      return result;
    }

看起來好像要執行更多東西，但是使用這個function時花費的時間瞬間下降到0.018S(4個thread)
看來果然還是要考慮的race condition尤其是在使用lib定義的function時
這裡看一下rand_r()的呼叫方式

unsigned int seed = time(NULL);
x = (double)rand_r(&seed)/(RAND_MAX + 1.0);
y = (double)rand_r(&seed)/(RAND_MAX + 1.0);

這裡宣告的東西都是thread local的，所以可以看成都是private的，只是對於rand_r()就算我們把seed變成global讓每個thread共用也沒關係，因為他們使用seed的時候只有"讀取"並沒有要改寫seed，他最後return的變數叫做result，所以沒有race condition要考慮

那local跟global的差異是?其實就是變數是否會繼承上一次的亂數結果而已

如果有兩個thread共用同一個seed產生rand_r亂數 2 ,3,5,7,11,13,17，那麼因為seed共用，兩個thread會各自分配得到這些亂數看OS給他們執行的順序決定

thread 1                thread 2
           <---  2
                 3 --->
           <---  5
                 7 --->
           <--- 11
                13 --->
           <--- 17

如果seed是private，那a、b代表兩個thread各自自己的seed

thread 1                thread 2
           <---  2a
                 2b --->
           <---  3a
                 3b --->
           <---  5a
                 5b --->
                 ::

那他們應該會得到一樣的結果不管OS給這些thread執行的順序是甚麼，不過對於我們這個作業來說，解決race condition就已經十分足夠了。