Параллельная быстрая сортировка на Хаскеле

от автора

Прим. перев.: Это перевод истории о том, как нелегко оказалось написать параллельную быструю сортировку (quicksort) на Хаскеле. Оригинал статьи написан в 2010 году, но, мне кажется, он до сих пор поучительный и во многом актуальный.

Есть много примеров того, как Хаскель делает простые проблемы сложными. Вероятно, самый известный из них—это решето Эратосфена, которое легко написать на любом императивном языке, но настолько сложно написать на Хаскеле, что почти все решения, которые преподавались в университетах и использовались в исследованиях последние 18 лет, оказались неправильными. На их несостоятельность обратила внимание Мелисса О’Нил [Melissa O’Neill] в своей важной научной работе "Настоящее решето Эратосфена". В ней приводится прекрасное описание того, что не так в старых подходах, и как их надо исправить. Решением Мелиссы было использовать очередь с приоритетом [priority queue] для реализации решета. Правильное решение оказалось в 10 раз длиннее, чем намного более простое решение на F# и в целых 100 раз длиннее, чем оригинальный изуродованный алгоритм на Хаскеле.

Сегодня быстрая сортировка—это новое решето Эратосфена. И программисты на Хаскеле снова обошли неспособность языка выразить этот алгоритм уродованием последнего. Новый вариант медленнее на порядки, но зато его можно легко записать на Хаскеле.

qsort []     = [] qsort (x:xs) = qsort (filter (< x) xs) ++ [x] ++ qsort (filter (>= x) xs) 

Этот код совершенно не соотносится с сущностью настоящего алгоритма быстрой сортировки, которая делает его таким эффективным (см. оригинальную статью Тони Хоара 1962 года о быстрой сортировке). А именно, перегруппировка [partitioning] массива без дополнительного выделения памяти [in-place partitioning using swaps].

Столкнувшись с проблемой написания паралелльной быстрой сортировки общего назначения на Хаскеле, Джим Эппл [Jim Apple] (который пишет кандидатскую по Хаскелю в Калифорнийском университете в Дейвисе, UC Davis) дал старт делу, написав следующий код:

import Data.HashTable as H import Data.Array.IO import Control.Parallel.Strategies import Control.Monad import System  exch a i r =     do tmpi <- readArray a i        tmpr <- readArray a r        writeArray a i tmpr        writeArray a i tmpi  bool a b c = if c then a else b  quicksort arr l r =   if r <= l then return () else do     i <- loop (l-1) r =<< readArray arr r     exch arr i r     withStrategy rpar $ quicksort arr l (i-1)     quicksort arr (i+1) r   where     loop i j v = do       (i', j') <- liftM2 (,) (find (>=v) (+1) (i+1)) (find (<=v) (subtract 1) (j-1))       if (i' < j') then exch arr i' j' >> loop i' j' v                    else return i'     find p f i = if i == l then return i                  else bool (return i) (find p f (f i)) . p =<< readArray arr i  main =      do [testSize] <- fmap (fmap read) getArgs        arr <- testPar testSize        ans <- readArray arr  (testSize `div` 2)        print ans  testPar testSize =     do x <- testArray testSize        quicksort x 0 (testSize - 1)        return x  testArray :: Int -> IO (IOArray Int Double) testArray testSize =      do ans <- newListArray (0,testSize-1) [fromIntegral $ H.hashString $ show i | i <- [1..testSize]]        return ans 

Этот алгоритм использует параллельные "стратегии" Хаскеля. Эта концепция была разработана, чтоб дать программистам на Хаскеле больше контроля над параллелизацией, но оказалось, что в единственной доступной имплементации течёт память и никому не удалось заставить её работать в этом коде: решение Джима содержит ошибку многопоточности [concurrency], из-за которой оно возвращает неправильные результаты почти при каждом вызове.

Тогда Пикер [Peaker] предложил свое решение на Хаскеле:

import Data.Array.IO import Control.Monad import Control.Concurrent  bool t _f True = t bool _t f False = f  swap arr i j = do   (iv, jv) <- liftM2 (,) (readArray arr i) (readArray arr j)   writeArray arr i jv   writeArray arr j iv  parallel fg bg = do   m <- newEmptyMVar   forkIO (bg >> putMVar m ())   fg >> takeMVar m  sort arr left right = when (left < right) $ do   pivot <- read right   loop pivot left (right - 1) (left - 1) right   where     read = readArray arr     sw = swap arr     find n pred i = bool (find n pred (n i)) (return i) . pred i =<< read i     move op d i pivot = bool (return op)                         (sw (d op) i >> return (d op)) =<<                         liftM (/=pivot) (read i)     loop pivot oi oj op oq = do       i <- find (+1) (const (<pivot)) oi       j <- find (subtract 1) (\idx cell -> cell>pivot && idx/=left) oj       if i < j         then do           sw i j           p <- move op (+1) i pivot           q <- move oq (subtract 1) j pivot           loop pivot (i + 1) (j - 1) p q         else do           sw i right           forM_ (zip [left..op-1] [i-1,i-2..]) $ uncurry sw           forM_ (zip [right-1,right-2..oq+1] [i+1..]) $ uncurry sw           let ni = if left >= op then i + 1 else right + i - oq               nj = if right-1 <= oq then i - 1 else left + i - op           let thresh = 1024               strat = if nj - left < thresh || right - ni < thresh                       then (>>)                       else parallel           sort arr left nj `strat` sort arr ni right  main = do   arr <- newListArray (0, 5) [3,1,7,2,4,8]   getElems arr >>= print   sort (arr :: IOArray Int Int) 0 5   getElems arr >>= print 

Это решение тоже оказалось с багами. Во-первых, оно содержит более тонкий баг многопоточности [concurrency], который приводит к неверным результатам лишь изредка. Пикер исправил этот баг в следующем коде:

import System.Time import System.Random import Data.Array.IO import Control.Monad import Control.Concurrent import Control.Exception import qualified Data.List as L  bool t _ True = t bool _ f False = f  swap arr i j = do   (iv, jv) <- liftM2 (,) (readArray arr i) (readArray arr j)   writeArray arr i jv   writeArray arr j iv  background task = do   m <- newEmptyMVar   forkIO (task >>= putMVar m)   return $ takeMVar m  parallel fg bg = do   wait <- background bg   fg >> wait  sort arr left right = when (left < right) $ do   pivot <- read right   loop pivot left (right - 1) (left - 1) right   where     read = readArray arr     sw = swap arr     find n pred i = bool (find n pred (n i)) (return i) . pred i =<< read i     move op d i pivot = bool (return op)                         (sw (d op) i >> return (d op)) =<<                         liftM (/=pivot) (read i)     swapRange px x nx y ny = if px x then sw x y >> swapRange px (nx x) nx (ny y) ny else return y     loop pivot oi oj op oq = do       i <- find (+1) (const (<pivot)) oi       j <- find (subtract 1) (\idx cell -> cell>pivot && idx/=left) oj       if i < j         then do           sw i j           p <- move op (+1) i pivot           q <- move oq (subtract 1) j pivot           loop pivot (i + 1) (j - 1) p q         else do           sw i right           nj <- swapRange (<op) left (+1) (i-1) (subtract 1)           ni <- swapRange (>oq) (right-1) (subtract 1) (i+1) (+1)           let thresh = 1024000               strat = if nj - left < thresh || right - ni < thresh                       then (>>)                       else parallel           sort arr left nj `strat` sort arr ni right  timed act = do   TOD beforeSec beforeUSec <- getClockTime   x <- act   TOD afterSec afterUSec <- getClockTime   return (fromIntegral (afterSec - beforeSec) +           fromIntegral (afterUSec - beforeUSec) / 1000000000000, x)  main = do   let n = 1000000   putStrLn "Making rands"   arr <- newListArray (0, n-1) =<< replicateM n (randomRIO (0, 1000000) >>= evaluate)   elems <- getElems arr   putStrLn "Now starting sort"   (timing, _) <- timed $ sort (arr :: IOArray Int Int) 0 (n-1)   print . (L.sort elems ==) =<< getElems arr   putStrLn $ "Sort took " ++ show timing ++ " seconds" 

Это решение действительно работает на маленьких входных массивах, но увеличение размера массива до 1 000 000 элементов приводит к переполнению стека. Было сделано две попытки проанализировать этот баг, здесь и здесь, но обе оказались неправильными. На самом деле, это баг в функции getElems стандартной библиотеки Хаскеля, которая переполняет стек на больших массивах.

Как ни странно, исправления еще нескольких багов, судя по всему, привели к реализации первой в мире параллельной быстрой сортировки общего назначения, написанной на Хаскеле. Более того, финальное решение на Хаскеле всего примерно на 55% медленнее, чем эквивалентное решение на F#. Будьте внимательны, это решение требует последнюю версию GHC, которая была выпущена несколько недель назад (прим. перев.: статья 2010 года, так что читателю беспокоиться не о чем).

Первые комментарии к оригинальной статье

Ganesh Sittampalam:

Поздравляю с обучением тому, как делать fork и synchronise в Хаскеле!

Jon Harrop (автор оригинала):

Поздравляю с проверкой твоей теории о том, что это будет «тривиально»…

ссылка на оригинал статьи https://habrahabr.ru/post/317348/


Комментарии

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *