M.Hiroi's Home Page
http://www.geocities.jp/m_hiroi/

Julia Language Programming

並行プログラミング超入門

[ Home | Light | Julia ]

はじめに

「並行 (concurrent) プログラミング」は複数のプログラムを並行に実行しますが、このとき複数の CPU で同時に動かす場合と、ひとつの CPU で複数のプログラムを動かす場合があります。一般的には、前者を「並列 (parallel) プログラミング」といい、複数のハードウェアを並列に実行することによる処理速度の向上が主な目的となります。

後者の場合、一定時間毎に実行するプログラムを切り替えることで、複数のプログラムを並列に実行しているかのように見せることができます。この処理を「時分割 (time sharing)」もしくは「タイム・スライス (time slice)」といいます。一般に、タイム・スライスは OS でサポートされている機能です。OS が実行するプログラムのことを「プロセス (process)」または「タスク (task)」といいます。Julia ではコルーチンのことをタスクというので、本稿ではプロセスと記述することにします。

並列的に動作するプログラムをひとつのプロセスだけで作るのはけっこう大変です。そこで、プロセス内では逐次的な処理にとどめ、複数のプロセス間で情報交換を行うことにより、全体で並列的な動作を実現することを考えます。このほうが自然にプログラムを記述できる場合があるのです。これが後者の主な目的となります。

プロセスは互いに独立したプログラムですが、OS にはプロセス間でデータをやり取りする機能 (プロセス間通信) が用意されています。たとえば、UNIX ライクな OS では「パイプ (pipe)」を使って複数のプログラム (コマンド) を連結することができます。この場合、パイプを通してデータがプログラムに送られ、各プログラムは並行に動作することになります。入出力処理で待たされるコマンドがあったとしても、その間に他のコマンドが実行されるので、各コマンドを逐次的に実行するよりも、効率的に処理することが可能です。

最近では、ひとつのプログラムの中で独立した複数の処理を記述できるようになりました。この機能を「スレッド (thread)」とか「マルチスレッド」いいます。スレッドは「縫い糸」という意味ですが、プログラムでは「制御の流れ」という意味で使われています。並列的な動作をプログラムする場合、逐次的な処理をひとつのスレッドに割り当て、複数のスレッドを並行に動作させることにより、全体で並列的な動作を実現するわけです。

一般的なスレッドは、一定時間毎に実行するスレッドを強制的に切り替えます。このとき、スレッドのスケジューリングは処理系が行います。これを「プリエンプティブ (preemptive)」といいます。コルーチンの場合、プログラムの実行は一定時間ごとに切り替わるものではなく、プログラム自身が実行を中断しないといけません。これを「ノンプリエンプティブ (nonpreemptive)」といいます。

また、スレッドは同じプロセス内に存在するので、メモリ空間を共有することができます。これを「共有メモリ」といいます。スレッド間の通信は共有メモリを使って簡単に行うことができますが、プリエンプティブなスレッドの場合、共有メモリのアクセス時に発生する「競合」が問題になります。このため、マルチスレッドをサポートしているプログラミング言語では、競合を回避するための仕組みが用意されています。

今回はコルーチンを使って擬似的なマルチプロセスを実現するプログラムを作ってみましょう。このプログラムはコルーチンを利用したノンプリエンプティブなものなので、競合について考慮する必要はありません。ただし、複数のプロセス間で通信を行うので、待ち合わせが必要になることがあります。この処理を「同期 (synchronization)」といいます。並行プログラミングの場合、通信と同期という 2 つの処理が重要になります。

●簡単なマルチプロセスの作成

それではプログラムを作りましょう。一般に、プロセスには優先順位 (プライオリティ) があるのですが、今回はプログラムを簡単にするため、すべてのプロセスは同じ優先順位とします。この場合、コルーチンをそのままプロセスとして扱うと簡単です。

最初に、メインプロセスをひとつ用意し、そこでコルーチンを生成して呼び出します。中断したコルーチンはキューに格納しておけばいいでしょう。つまり、キューからコルーチンを取り出して実行を再開し、中断したら再びキューに追加すればいいわけです。コルーチンの実行が終了した場合、そのコルーチンはキューに追加しません。これで生成したプロセスを擬似的にですが並行に実行することができます。

プログラムは次のようになります。

リスト : 簡単なマルチプロセス

# キュー
const queue = Task[]

# 挿入
enqueue(p) = push!(queue, p)

# 削除
dequeue() = shift!(queue)

# キューは空か
is_empty() = queue == []

# プロセスの生成
fork(fn) = enqueue(Task(fn))

# メインプロセス
function main_process(args...)
    # プロセスの登録
    for fn in args; fork(fn); end
    # 実行
    while !is_empty()
        p = dequeue()
        if consume(p); enqueue(p); end
    end
end

# 実行権を放棄するだけ
yield() = produce(true)

大域変数 queue には中断したプロセスを格納する配列をセットします。これをキューとして使います。プロセスの生成は関数 fork() で行います。引数 fn は引数なしの関数とします。プロセスの終了は false で表すことにします。あとは Task() でコルーチンを生成し、それを関数 enqueue() でキューに追加します。

メインプロセスの実行は関数 main_process() で行います。引数はプロセスの実体を表す関数です。最初に引数 args から関数をひとつずつ取り出し、それを fork() に渡してプロセスを生成します。あとはキューからプロセスを順番に取り出して変数 p にセットし、consume() でプロセス p の実行を再開します。返り値が真の場合、プロセス p はまだ終了していないので、p をキューに追加します。返り値が偽の場合はキューに追加しません。

関数 yield() はプロセスの実行権を他のプロセスに渡すだけです。これは produce() で true を返すだけです。これで main_process() に戻り、次のプロセスが実行されます。

●簡単な実行例

それでは簡単な例を示しましょう。次のリストを見てください。

リスト : 簡単なテスト

function test0(name, n)
    while n > 0
        println("$name $n")
        n -= 1
        yield()
    end
    false
end

main_process(() -> test0("foo", 5),
             () -> test0("bar", 6),
             () -> test0("baz", 4))

関数 test0() は name を n 回表示します。name と n を表示したあと、yield() で実行権を放棄します。これで他のプロセスが実行されるので、複数のプロセスを並行に動作させることができます。実行例を示します。

foo 5
bar 6
baz 4
foo 4
bar 5
baz 3
foo 3
bar 4
baz 2
foo 2
bar 3
baz 1
foo 1
bar 2
bar 1

このように、他のプロセスと通信を行わない場合、プログラムはとても簡単になります。

●schedule() の使い方

実をいうと、Julia には main_process() に相当する関数 schedule() またはマクロ @shcedule() が用意されています。schedule() はタスクをスケジューラのキューに追加して、システムがアイドルなときにタスクを実行します。スケジューラーはメインのプログラムと並行に動作していると考えてください。このため、メインのプログラムが終了すると、スケジューラで実行されているタスクも終了します。ご注意ください。

簡単な実行例を示します。

リスト : schedule() の使い方

function test1(name, n)
    while n > 0
        println("$name $n")
        n -= 1
        yield()
    end
    false
end

t1 = @task test1("foo", 5)
t2 = @task test1("bar", 6)
t3 = @task test1("baz", 4)
schedule(t1)
schedule(t2)
schedule(t3)
wait(t2)

test1() の yield() は Julia の関数で、タスクの実行を中断して制御をスケジューラに戻します。スケジューラのキューに他のタスクがあれば、それをキューから取り出して実行します。一番最後の wait() は、引数がタスクの場合、それが終了するまでプログラムの実行を一時停止します。これがないとメインのプログラムが先に終了するので、スケジューラに登録したタスクは実行されません。

実行結果は test0() と同じです。

foo 5
bar 6
baz 4
foo 4
bar 5
baz 3
foo 3
bar 4
baz 2
foo 2
bar 3
baz 1
foo 1
bar 2
bar 1

●チャネルによる同期処理

次はプロセス間で通信を行う場合を考えてみましょう。この場合、いろいろな方法が考えられますが、今回は Julia に用意されている「Channel (チャネル)」というデータ型を使ってみましょう。チャネルの生成はコンストラクタ Channel() を使います。

Channel{T}(size::Int)  # T 型のチャネル (大きさ size)
Channel(size::Int)     # Any 型のチャネル (大きさ size)
Channel()              # Any 型のチャネル (大きさ 32)

チャネルの基本はキューと同じですが、データを追加する関数 put!() とデータを取り出す関数 take!() で待ち合わせを行うところが異なります。put!() では、キューが満杯のときに待ち合わせを行います。逆に take!() の場合、キューが空のときに待ち合わせを行います。キューの大きさが少ない場合でも、データを書き込むプロセスと読み出すプロセスが並行に動作することで、キューの大きさ以上のデータを受け渡すことができます。

それでは簡単な実行例を示します。次のリストを見てください。

リスト : チャネルの実行例

# チャネル
const ch = Channel(4)

# データを送る
function send_color(color, n)
    for i in 1 : n
        put!(ch, color)
    end
end

# データを受け取る
function receive_color(n)
    for i in 1 : n
        println(take!(ch))
    end
end

# テスト
@schedule send_color("red", 8)
@schedule send_color("blue", 8)
@schedule send_color("yellow", 8)
t1 = @task receive_color(24)
schedule(t1)
wait(t1)

Channel(4) でチャネルを生成して大域変数 ch に格納します。大きさは 4 とします。send_color() はデータ (color) を n 個チャネルに書き込みます。receive_color() は n 個のデータをチャネルから取り出して表示します。テストは、キューに書き込むプロセスを 3 つ生成し、取り出すプロセスを 1 つ生成します。データを書き込む総数は 8 * 3 = 24 個なので、取り出すデータ数も 24 個とします。

それでは実行結果を示します。

red
red
red
red
red
blue
yellow
red
blue
yellow
red
blue
yellow
red
blue
yellow
blue
yellow
blue
yellow
blue
yellow
blue
yellow

24 個のデータすべて表示されていますね。

●哲学者の食事

最後に、「哲学者の食事」という並行プログラミングでは有名な問題を解いてみましょう。

[哲学者の食事]

5 人の哲学者が丸いテーブルに座っています.テーブルの中央にはスパゲッティが盛られた大皿があり、哲学者の間には 5 本のフォークが置かれています。哲学者は思索することとスパゲッティを食べることを繰り返します。食事のときには 2 本のフォークを持たなければなりません。食事が終わると 2 本のフォークを元の位置に戻します。

詳しい説明は 食事する哲学者の問題 -- Wikipedia をお読みください。

それではプログラムを作りましょう。最初にフォークを操作する関数を定義します。

リスト : フォークを操作する関数

# フォーク
const forks = [true, true, true, true, true]

# フォークの番号を求める
function fork_num(person, side)
    n = person
    if side == :left
        n = n + 1
        if n > length(forks)
            n = 1
        end
    end
    n
end

# フォークがあるか?
function is_fork(person, side)
    forks[fork_num(person, side)]
end

# フォークの書き換え
function fork_set(person, side, val)
    forks[fork_num(person, side)] = val
end

# フォークを取る
function get_fork(person, side)
    # 待ち合わせ
    while true
        if is_fork(person, side); break; end
        yield()
    end
    fork_set(person, side, false)
    yield()
end

# フォークを置く
function put_fork(person, side)
    fork_set(person, side, true)
    yield()
end

フォークの有無は真偽値で表して、それを配列に格納します。配列は大域変数 forks にセットします。n 番目の哲学者の場合、右側のフォークがベクタの n 番目の要素、左側のフォークが n + 1 番目の要素になります。関数 fork_num() はフォークに対応する番号を返します。関数 is_fork() は n 番目の哲学者の side にフォークがあるとき真を返します。fork_set() は n 番目の哲学者の side 側にあるフォークの値を val で書き換えます。

get_fork() はフォークを取る関数です。while ループで is_fork() が真を返すまで待ちます。このとき、yield() を呼び出して他のプロセスに実行権を渡すことをお忘れなく。そのあとで、fork_set() で対応するフォークの値を false() に書き換えます。put_fork() はフォークを元に戻す関数です。fork_set() でフォークの値を true に書き換えます。

今回はノンプリエンプティブなコルーチンを使用しているので、forks をアクセスするときに競合は発生しません。プリエンプティブなマルチスレッドを使用する場合、共有メモリにアクセスするときは競合について考慮する必要があります。ご注意ください。

次は哲学者の動作をプログラムします。次のリストを見てください。

リスト : 哲学者の動作

# 哲学者の動作
function person0(n)
    m = 0
    while m < 2
        println("Philosopher $n is thinking")
        sleep(1)
        get_fork(n, :right)
        get_fork(n, :left)
        println("Philosopher $n is eating")
        sleep(1)
        put_fork(n, :right)
        put_fork(n, :left)
        m = m + 1
    end
    println("Philosopher $n is sleeping")
end

関数 person0 の引数 n は哲学者の番号を表します。変数 m は食事をした回数です。2 回食事をしたら処理を終了します。食事をする場合、最初に get_fork() で右側のフォークを取り、次に左側のフォークを取ります。食事を終えたら put_fork() で右側のフォークを置き、次に左側のフォークを置きます。

thinking と eating のあとは、関数 sleep() でプログラムの実行を休止します。引数の単位は秒 (sec) です。ミリ秒は浮動小数点数で指定してください。なお、スケジューラに登録されたプログラムで sleep() を使用すると、実行権は他のプログラムに移ります。ご注意くださいませ。

このように、マルチプロセスを使うと簡単にプログラムできますが、実は並行プログラミング特有の大きな問題点があるのです。これはプログラムを実行してみるとわかります。

●実行結果 (1)

プログラムの実行は関数 test() で行います。

リスト : 実行

function test(person)
    ps = map(x -> Task(() -> person(x)), [1, 2, 3, 4, 5])
    for p in ps; schedule(p); end
    for p in ps; wait(p); end
end

test(person0)

map() でプロセスを生成し、それを schedule() でスケジューラに登録します。あとは、wait() でプロセスが終了するのを待つだけです。実行結果は次のようになります。

Philosopher 1 is thinking
Philosopher 2 is thinking
Philosopher 3 is thinking
Philosopher 4 is thinking
Philosopher 5 is thinking

このように、すべてのプロセスが待ち状態となり進むことができなくなります。これを「デッドロック (deadlock)」といいます。この場合、1 番目の哲学者の右側のフォークは、5 番目の哲学者の左側のフォークになります。各哲学者が右側のフォークを取り、左側のフォークが置かれるのを待つときにデッドロックとなるわけです。

●デッドロックの防止

デッドロックを防止する簡単な方法は、右側のフォークを取っても左側のフォークを取れないときは、右側のフォークを元に戻すことです。プログラムは次のようになります。

リスト : デッドロックの防止 (1)

function person1(n)
    m = 0
    while m < 2
        println("Philosopher $n is thinking")
        sleep(1)
        get_fork(n, :right)
        if is_fork(n, :left)
            get_fork(n, :left)
            println("Philosopher $n is eating")
            sleep(1)
            put_fork(n, :right)
            put_fork(n, :left)
            m = m + 1
        else
            put_fork(n, :right)
        end
    end
    println("Philosopher $n is sleeping")
end

右側のフォークを取ったあと、is_fork() で左側のフォークをチェックします。フォークがある場合は、fork_set() で値を true に書き換えます。これでフォークを取って食事をすることができます。左側のフォークがない場合は put_fork() で右側のフォークを元に戻します。

●実行結果 (2)

実行結果は次のようになります。

Philosopher 1 is thinking
Philosopher 2 is thinking
Philosopher 3 is thinking
Philosopher 4 is thinking
Philosopher 5 is thinking
Philosopher 1 is thinking
Philosopher 2 is thinking
Philosopher 3 is thinking
Philosopher 4 is thinking
Philosopher 5 is eating
Philosopher 2 is thinking
Philosopher 3 is thinking
Philosopher 4 is eating
Philosopher 5 is thinking
Philosopher 1 is eating
Philosopher 3 is eating
Philosopher 4 is thinking
Philosopher 1 is thinking
Philosopher 5 is eating
Philosopher 2 is thinking
Philosopher 3 is thinking
Philosopher 1 is thinking
Philosopher 5 is sleeping
Philosopher 2 is eating
Philosopher 4 is eating
Philosopher 1 is thinking
Philosopher 2 is thinking
Philosopher 4 is sleeping
Philosopher 3 is eating
Philosopher 1 is eating
Philosopher 3 is sleeping
Philosopher 1 is sleeping
Philosopher 2 is eating
Philosopher 2 is sleeping

このように、デッドロックしないで実行することができます。

●デッドロックの防止 (2)

もうひとつ簡単な方法を紹介しましょう。奇数番目の哲学者は、まず左側のフォークを取り上げてから右側のフォークを取り、偶数番目の哲学者は、今までのように右側のフォークを取り上げてから左側のフォークを取ります。こんな簡単な方法で動作するのは不思議なように思います。たとえば、哲学者が 2 人の場合を考えてみてください。

哲学者 0 の右側のフォークを A、左側のフォークを B とします。哲学者 1 からみると、B が右側のフォークで、A が左側のフォークになります。デッドロックは、哲学者 0 が A を取り、哲学者 1 が B を取ったときに発生します。ここで、哲学者 1 が左側のフォーク A から取るようにします。先に哲学者 0 が A を取った場合、哲学者 1 は A があくまで待つことになるので、哲学者 0 はフォーク B を取って食事をすることができます。哲学者 1 が先にフォーク A を取った場合も同じです。これでデッドロックを防止することができます。

プログラムは次のようになります。

リスト : デッドロックの防止 (2)

function person2(n)
    m = 0
    while m < 2
        println("Philosopher $n is thinking")
        sleep(1)
        if n % 2 == 0
            get_fork(n, :right)
            get_fork(n, :left)
        else
            get_fork(n, :left)
            get_fork(n, :right)
        end
        println("Philosopher $n is eating")
        sleep(1)
        put_fork(n, :right)
        put_fork(n, :left)
        m = m + 1
    end
    println("Philosopher $n is sleeping")
end

if 文で n が偶数の場合は右側から、奇数の場合は左側のフォークから取るように処理を分けるだけです。

●実行結果 (3)

実行結果は次のようになります。

Philosopher 1 is thinking
Philosopher 2 is thinking
Philosopher 3 is thinking
Philosopher 4 is thinking
Philosopher 5 is thinking
Philosopher 3 is eating
Philosopher 5 is eating
Philosopher 3 is thinking
Philosopher 5 is thinking
Philosopher 1 is eating
Philosopher 4 is eating
Philosopher 1 is thinking
Philosopher 4 is thinking
Philosopher 3 is eating
Philosopher 5 is eating
Philosopher 2 is eating
Philosopher 4 is eating
Philosopher 3 is sleeping
Philosopher 5 is sleeping
Philosopher 2 is thinking
Philosopher 4 is sleeping
Philosopher 1 is eating
Philosopher 1 is sleeping
Philosopher 2 is eating
Philosopher 2 is sleeping

正常に動作していますね。興味のある方はいろいろ試してみてください。

●並列プログラミング

今度は Julia で「並列 (parallel) プログラミング」に挑戦してみましょう。なお、M.Hiroi は Julia もそうですが並列プログラミングに関しても初心者にすぎません。なにぶんにも初心者が作るプログラムなので、何か問題点や間違いがあるかもしれません。お気づきの点がありましたら、ご指摘いただけると助かります。

なお、並列プログラミングの基本は拙作のページ Julia の基礎知識: 並列プログラミング で簡単に説明していますので、そちらをお読みくださいませ。

●数値積分

まず最初に、数値積分で円周率πを求めてみましょう。区間 [a, b] の定積分∫f(x)dx を数値的に求めるには、区間を細分して小区間の面積を求めて足し上げます。小区間の面積を求める一番簡単な方法は長方形や台形で近似することです。

  1. (b - a) * f(a)
  2. (b - a) * f(b)
  3. (b - a) * f((a + b) / 2)
  4. (1/2) * (b - a) * (f(a) + f(b))

1 は左端の値 f(a) を、2 は右端の値 f(b) を、3 は中間点の値 f((a + b) / 2) を使って長方形の面積を計算します。1, 2, 3 の中では 3 番目の方法が一番精度が高く、これを「中点則」といいます。4 は台形で面積を計算するので「台形則」といいます。中点則は拙作のページ 数値積分 でプログラムしたことがあるので、今回は台形則を使ってみましょう。

πは次の式で求めることができます。

      1
π = ∫(4 / (1 + x * x)) dx
      0

プログラムは次のようになります。

リスト : 台形則

function trapezoid(n)
    const w = 1.0 / n
    f(x) = 4.0 / (1.0 + x * x)
    s = 0.0
    for i in 0 : n - 1
        s += 0.5 * w * (f(i * w) + f((i + 1) * w))
    end
    s
end

for n in 1 : 8
    x = 10 ^ n
    p = trapezoid(x)
    e = pi - p
    println("--- $x ---")
    println("$p $e")
end
--- 10 ---
3.1399259889071587 0.0016666646826344333
--- 100 ---
3.141575986923131 1.666666666233496e-5
--- 1000 ---
3.1415924869231278 1.6666666535769536e-7
--- 10000 ---
3.1415926519231183 1.6666747981730623e-9
--- 100000 ---
3.141592653573132 1.66608948859448e-11
--- 1000000 ---
3.1415926535897176 7.549516567451064e-14
--- 10000000 ---
3.1415926535904326 -6.394884621840902e-13
--- 100000000 ---
3.1415926535893615 4.3165471197426086e-13

台形則や中点則の場合、分割数を 10 倍すると誤差はほぼ 1/100 になります。ただし、浮動小数点数 (Float64) の計算には誤差があるので、精度には限界があります。台形則の場合、分割数を 1,000,000 より増やしても精度は高くなりませんが、並列プログラミングの例題として、あえて分割数を 10,000,000 に設定して試してみましょう。

プログラムは次のようになります。

リスト : 台形則の並列化

@everywhere function trapezoid_para(n, d)
    const w = 1.0 / n
    f(x) = 4.0 / (1.0 + x * x)
    function trapezoid(n, m)
        s = 0.0
        for i in n : m - 1
            s += 0.5 * w * (f(i * w) + f((i + 1) * w))
        end
        s
    end
    k = div(n, d)
    @parallel (+) for m in 0 : k : n - 1
        trapezoid(m, m + k)
    end
end

M.Hiroi のパソコンは物理コア数が 4 で、1 コアにつきハイパースレッディングで 2 分割できるので、CPU_CORES の値は 4 * 2 = 8 になります。分割数を 1, 2, 4, 8 として、実行時間を計測しました。なお、実行時間は 5 回計測して、最大値と最小値を除いた 3 つの値の平均値を計算したものです。

  表 : 実行結果

 d : 時間 : 効率
---+------+------
 1 : 2.33 : 1.00
 2 : 1.17 : 1.99
 4 : 0.78 : 2.98
 8 : 0.70 : 3.32

実行環境 : Windows 7, Core i7-2670QM 2.20GHz

Go 言語入門 並列プログラミング で同じようなプログラムを作ったことがありますが、Julia の場合、物理コア数を超えると並列化の効果は少なくなるようです。今後のバージョンアップで並列処理のパフォーマンスが改善されることを期待したいと思います。

●モンテカルロ法

乱数を使って数学や物理などの問題を解くアルゴリズムを「モンテカルロ法 (Monte Carlo methods)」といいます。簡単な例題として、円周率πをモンテカルロ法で求めてみましょう。

正方形の領域 (0 <= x < 1, 0 <= y < 1) に乱数で点を打ちます。乱数であれば点は領域内に一様に分布するので、x2 + y2 < 1 の円内に入る確率は π / 4 になります。つまり、(円内の点の個数 / 点の総数) の値は 0.7853... になるはずです。たくさん点を打つほど値は π / 4 に近づくはずですが、コンピュータの乱数は疑似乱数なので規則性が生じてしまい、値の精度にはどうしても限界があります。

また、たくさん点を打つと実行時間がかかるようになりますが、並列に処理することで実行時間を短縮することができます。たとえば、100 万個の点を打ってπを求める処理を並列に 8 回行って平均値を計算すれば、800 万個の点を打ってπを求めたことと同じになります。

●乱数列の生成

Julia の場合、関数 rand() で乱数を取得することができます。Julia の rand() は高機能です。

rand([rng][, S][, dims...])

rng は乱数発生ジェネレータです。省略すると大域変数に格納されたジェネレータが使用されます。S はデータ型やコレクション (['a', 'b', 'c'] とか 1 : 10 など) を指定します。データ型 S が整数の場合、乱数の範囲は [typemin(S), typemax(S)] になります。浮動小数点数の場合は [0, 1.0) (0 以上 1.0 未満) の範囲になります。dims を指定すると、dims で指定した大きさの配列に乱数を格納して返します。

簡単な例を示しましょう。

julia> for _ in 1 : 10; println(rand()); end
0.5939723464136937
0.024220850690386264
0.648525189990292
0.22083957896326245
0.1309902585960312
0.47506648591383205
0.41566489255438244
0.3515387127792815
0.9575547284952766
0.7399959494399366

julia> rand(3,3)
3x3 Array{Float64,2}:
 0.184927  0.541715  0.957923
 0.169656  0.374786  0.588924
 0.875136  0.777696  0.418522

julia> for _ in 1 : 10; println(rand(1:100)); end
77
5
29
22
70
57
44
2
54
15

julia> rand(1:100, 4,4)
4x4 Array{Int32,2}:
 39  65  68  99
 77  98  75  97
 67  39  33  64
 58  75  87  26

乱数発生ジェネレータは MersenneTwister() で生成することができます。乱数を生成するアルゴリズムは名前が表すようにメルセンヌ・ツイスターが使われています。

MersenneTwister([seed])

一般に、乱数の元になるデータを「シード (seed : 種)」といいます。乱数ジェネレータを生成する seed の値を変えることで、異なる乱数列を発生させることができます。簡単な例を示しましょう。

julia> r0 = MersenneTwister(0)
・・・省略・・・

julia> rand(r0, 10)
10-element Array{Float64,1}:
 0.823648
 0.910357
 0.164566
 0.177329
 0.27888
 0.203477
 0.0423017
 0.0682693
 0.361828
 0.973216

julia> r1 = MersenneTwister(1)
・・・省略・・・

julia> rand(r1, 10)
10-element Array{Float64,1}:
 0.236033
 0.346517
 0.312707
 0.00790928
 0.488613
 0.210968
 0.951916
 0.999905
 0.251662
 0.986666

julia> r = MersenneTwister()
・・・省略・・・

julia> rand(r, 10)
10-element Array{Float64,1}:
 0.823648
 0.910357
 0.164566
 0.177329
 0.27888
 0.203477
 0.0423017
 0.0682693
 0.361828
 0.973216

seed を省略すると 0 を指定したのと同じ乱数列が生成されます。

●モンテカルロ法の並列化

Julia の場合、大域変数は各ワーカープロセスごとに用意されます。大域的な乱数ジェネレータも同じです。シードを関数 srand(seed) で設定すれば、ワーカープロセスごとに異なる乱数列を発生することができますが、今回は乱数の例題もかねて、プロセスごとにジェネレータを生成することにします。

プログラムは次のようになります。

リスト : モンテカルロ法

@everywhere function monte_pi(n, s)
    rng = MersenneTwister(s)       # 大域的な乱数ジェネレータを使用する場合は
    c = 0                          # srand(s) でシードを設定して、
    for _ in 1 : n                 # x = rand()
        x = rand(rng)              # y = rand()
        y = rand(rng)              # として乱数を取得すればよい
        if x * x + y * y < 1.0
            c += 1
        end
    end
    4.0 * c / n
end

function test_monte_pi(n)
    for i in [1, 2, 4, 8]
        println("---- $i ----")
        m = div(n, i)
        p = @time @parallel (+) for s in 1 : i
                monte_pi(m, s)
            end
        println(p / i)
    end
end
julia> test_monte_pi(100000000)
---- 1 ----
  1.861125 seconds (1.01 k allocations: 58.609 KB)
3.14161456
---- 2 ----
  0.958163 seconds (1.97 k allocations: 115.629 KB)
3.1416066799999998
---- 4 ----
  0.512461 seconds (3.94 k allocations: 232.980 KB)
3.14147464
---- 8 ----
  0.438859 seconds (7.85 k allocations: 460.918 KB, 0.81% gc time)
3.14177268

実行環境 : Windows 7, Core i7-2670QM 2.20GHz

4 分割で 3.6 倍の高速化を達成できました。やっぱり物理コア数を超えると並列化の効果は落ちるようです。興味のある方は試してみてください。

●参考文献, URL

  1. Paul Graham (著),野田 開 (訳), 『On Lisp』, Web 版
  2. Timothy Buddy (著), 吉田雄二 (監修), 長谷川明生・大田義勝 (訳), 『Little Smalltake 入門』, アスキー出版, 1989
  3. Ravi Sethi (著), 神林靖 (訳), 『プログラミング言語の概念と構造』, アジソンウェスレイ, 1995

●プログラムリスト

リスト : 哲学者の食事

# フォーク
const forks = [true, true, true, true, true]

# フォークの番号を求める
function fork_num(person, side)
    n = person
    if side == :left
        n = n + 1
        if n > length(forks)
            n = 1
        end
    end
    n
end

# フォークがあるか?
function is_fork(person, side)
    forks[fork_num(person, side)]
end

# フォークの書き換え
function fork_set(person, side, val)
    forks[fork_num(person, side)] = val
end

# フォークを取る
function get_fork(person, side)
    # 待ち合わせ
    while true
        if is_fork(person, side); break; end
        yield()
    end
    fork_set(person, side, false)
    yield()
end

# フォークを置く
function put_fork(person, side)
    fork_set(person, side, true)
    yield()
end

# 哲学者の動作
function person0(n)
    m = 0
    while m < 2
        println("Philosopher $n is thinking")
        sleep(1)
        get_fork(n, :right)
        get_fork(n, :left)
        println("Philosopher $n is eating")
        sleep(1)
        put_fork(n, :right)
        put_fork(n, :left)
        m = m + 1
    end
    println("Philosopher $n is sleeping")
end

function person1(n)
    m = 0
    while m < 2
        println("Philosopher $n is thinking")
        sleep(1)
        get_fork(n, :right)
        if is_fork(n, :left)
            get_fork(n, :left)
            println("Philosopher $n is eating")
            sleep(1)
            put_fork(n, :right)
            put_fork(n, :left)
            m = m + 1
        else
            put_fork(n, :right)
        end
    end
    println("Philosopher $n is sleeping")
end

function person2(n)
    m = 0
    while m < 2
        println("Philosopher $n is thinking")
        sleep(1)
        if n % 2 == 0
            get_fork(n, :right)
            get_fork(n, :left)
        else
            get_fork(n, :left)
            get_fork(n, :right)
        end
        println("Philosopher $n is eating")
        sleep(1)
        put_fork(n, :right)
        put_fork(n, :left)
        m = m + 1
    end
    println("Philosopher $n is sleeping")
end

function test(person)
    ps = map(x -> Task(() -> person(x)), [1, 2, 3, 4, 5])
    for p in ps; schedule(p); end
    for p in ps; wait(p); end
end

#test(person0)
#test(person1)
test(person2)

Copyright (C) 2016 Makoto Hiroi
All rights reserved.

[ Home | Light | Julia ]