用 OpenMP 實作平行化的 Quicksort

OpenMP 和 Pthread 不同，我們沒有辦法很精確地控制有多少的 Thread 正在執行。所以如果要實作 Quicksort，我們不能直接使用 Recursion。相反地，我們必須引入 Threading Pool 的概念，讓各個 Thread 主動去某個工作清單拿工作。概念程式碼如下：

void quicksort_impl(sort_value_t array[],
                    stack_t *s,
                    size_t *busy_threads) {
  int idle = 1;

  size_t begin = 0;
  size_t end = 0;

  /* We should loop forever until everyone is idle */
  while (1) {
    while (idle) {
      /* This thread is idling.
         Try to get a new task from the stack. */

#pragma omp critical
      if (!stack_empty(s)) {
        /* Get the pending sort range, start work now! */
        sort_range_t *range = stack_top(s);
        begin = range->begin;
        end = range->end;
        stack_pop(s);

        ++*busy_threads;
        idle = 0;
      }

      /* If everyone finished his work, then leave now ... */
      if (*busy_threads == 0) {
        return;
      }
    }

    size_t size = end - begin;

    if (size <= 50) {
      sort_value_t *seg = array + begin;

      /* Use backward insertion sort when the sort range is small */
      size_t i, j;
      for (i = 1; i < size; ++i) {
        sort_value_t cur = seg[i];
        /* Move the number backward if they are greater than cur */
        for (j = i; j > 0 && seg[j - 1] > cur; --j) {
          seg[j] = seg[j - 1];
        }
        seg[j] = cur;
      }

      /* Mark begin as end so that we can acquire next sort range */
      begin = end;

#pragma omp critical
      --*busy_threads;

      idle = 1;

      /* Finished our range continue to acquire next range */
      continue;
    }


    /* Partition our range */
    sort_value_t *seg = array + begin;

    size_t pivot_index = rand() % size;
    SWAP(seg[0], seg[pivot_index]);
    sort_value_t pivot = seg[0];

    size_t left = 1;
    size_t right = size - 1;

    while (left < right) {
      while (seg[left] <= pivot && left < right) { left++; }
      while (seg[right] >= pivot && left < right) { right--; }
      SWAP(seg[left], seg[right]);
    }

    if (seg[left] > pivot) {
      left--;
    }
    SWAP(seg[0], seg[left]);

    /* Push left subrange to stack */
#pragma omp critical
    {
      sort_range_t range;
      range.begin = begin;
      range.end = begin + left;
      stack_push(s, &range);
    }

    /* Continue to sort right subrange */
    begin += left + 1;
  }
}

void quicksort(sort_value_t array[], size_t size) {
  size_t busy_threads = 0;
  stack_t *s = stack_new();
  sort_range_t all = { 0 , size };
  stack_push(s, &all);

  /* Run quicksort_impl in parallel */
#pragma omp parallel
  quicksort_impl(array, s, &busy_threads);

  stack_delete(s);
}

Haskell 與 State Monad

昨天試著用 Haskell 寫 State Monad，花了一點時間才寫出一個簡單的猜數字遊戲。感覺要在 Functional Programming Language 寫出 State 的概念比較困難...

import Control.Monad.State
import System.IO
import System.Random

gameLB, gameUB :: Integer
gameLB = 0
gameUB = 100

data GameState = GameState { steps :: Integer ,
                             lb :: Integer,
                             ub :: Integer }

initialState :: Integer -> Integer -> GameState
initialState l u = GameState { steps = 0, lb = l, ub = u }

updateState :: Integer -> Integer -> GameState -> GameState
updateState ans guess (GameState { steps = s, lb = l, ub = u }) =
    case (compare ans guess) of
      EQ -> GameState { steps = s + 1, lb = ans, ub = ans }
      LT -> GameState { steps = s + 1, lb = l, ub = min (guess - 1) u }
      GT -> GameState { steps = s + 1, lb = max l (guess + 1), ub = u }

range :: GameState -> (Integer, Integer)
range (GameState { lb = l, ub = u }) = (l , u)

stepCount :: GameState -> Integer
stepCount (GameState { steps = s }) = s

gameSession :: Integer -> StateT GameState IO Integer
gameSession ans =
    do (l, u) <- gets range
       lift $ putStr $ "Guess number (" ++ show l ++ ".." ++ show u ++ "): "
       lift $ hFlush stdout
       str <- lift $ getLine
       let num = read str
       modify (updateState ans num)
       case (compare ans num) of
         EQ -> gets stepCount >>= return
         LT -> gameSession ans >>= return
         GT -> gameSession ans >>= return

main :: IO ()
main = do rand <- getStdRandom (randomR (gameLB, gameUB))
          s <- evalStateT (gameSession rand) (initialState gameLB gameUB)
          putStrLn $ "Steps used: " ++ show s

備註：上面的程式用 ghc --make Guess.hs 就可以編譯了（有用到 mtl package）