並列計算ツールキット
The Parallel Computing Toolkit

異機種間ネットワーク上でMathematica を使って問題を並列処理
Solving problems with Mathematica in parallel on a heterogeneous network

Roman E. Maeder

MathConsult Dr. R. Mäder
Samstagernstraße 58a
8832 Wollerau
Switzerland
http://www.mathconsult.ch

MathLink 通信プロトコルは，複数のMathematica カーネルプロセスをMathematica 内部から制御するのに用いることができます．この機能により，並列プログラミングのための分散メモリ環境の実装が可能になります．並列計算用のMap等の並列言語用の命令は，これらの基本的な機能上に簡単に実現できます．

The MathLink communication protocol can be used to control several Mathematica kernel processes from within Mathematica. This feature allows the implementation of a distributed-memory environment for parallel programming. Parallel language constructs, such as a parallel version of Map, can easily be implemented on top of these primitive operations.

ここでは，Mathematica における並列プログラミングのライブラリについて述べ，仮想共有メモリや同期方法を含む，いくつかの使用例を示します．

We describe a library for parallel programming in Mathematica and give several examples of its use, including virtual shared memory and synchronization methods.

ライブラリは完全にMathematica で書かれているため，マシン非依存です．また，共通のファイルシステムを使わずに，異機種間ネットワークでも使うことができます．すべてのライブラリとアプリケーションコードは，MathLink を通して分散されます．

The library is written completely in Mathematica and is, therefore, machine-independent, allowing its use in heterogeneous networks and without common file systems. All library and application code is distributed through MathLink.

© 2001 MathConsult Dr. R. Mäder. All rights reserved.

はじめに
Introduction

並列計算（システムプログラマの視点）
Parallel computing ( a system programmer's view)

このコンテキストでは
In this context

並列計算ツールキットの機能
Features of the Parallel Computing Toolkit

MathLink

カーネルの制御
Controlling Kernels

ここでは初期設定を評価しないでください．

Don't evaluate initialization at this point

別のローカルカーネルを起動してください（これを行うコマンドは，オペレーティングシステムによって異なることがあります）．

Start another local kernel (the command to do so may depend on the operating system):

local = LinkLaunch["math -noinit -mathlink"]

rshかsshでリモートカーネルを起動してください．

Start a remote kernel with rsh or ssh:

remote=LinkLaunch["ssh1 vger.wolfram.com /home/rd2/scripts/math30 -mathlink -noinit"]

math -mathlink -linkcreateを選択した後，結果として表示されるポートにLinkConnectで接続することで，リモートカーネルを手動で開始することもできます．

You can also start the remote kernel manually with math -mathlink -linkcreate and then connect to the resulting port using LinkConnect.

初期行番号パケットを取り出します．

Swallowing initial line-number packet.

LinkRead[local]

LinkRead[remote]

ブロッキングを防ぎます．

Avoid blocking:

LinkReadyQ[local]

評価する式を送信します．

Sending expressions for evaluation:

LinkWrite[local, Unevaluated[EvaluatePacket[1 + 2]]]

結果をポーリングします．

Polling for result

If[LinkReadyQ[local], LinkRead[local], Null]

他のパケットと共に処理します．

Dealing with other packets

LinkWrite[remote, Unevaluated[EvaluatePacket[1/0]]]

If[LinkReadyQ[remote], LinkRead[remote], Null]

セッションを終了します．

Ending session

LinkClose /@ {local, remote}

配置ルーチン
Lining up Workers

MathLink との接続の開始・終了の高レベルルーチン（OSに依存）です．

Higher-level routines for starting/closing MathLink connections (OS dependent).

LaunchSlave["host", "cmdTemplate", options…]

ConnectSlave["port@host", options…]

cmdTemplate はStringFormで用いられる文字列で，`1`がホスト名で置き換えられます．2番目の形式は待機中のカーネルに接続します．この形式は，セキュリティのために直接リモートコマンドを実行することができないときに用いることができます．

cmdTemplate is a string used in StringForm, `1` will be replaced by the hostname. The second form connects to a waiting kernel; this form can be used if security prevents you from executing a remote command directly.

ローカルカーネルを起動します．

Start local kernel:

LaunchSlave["localhost", "math -mathlink"]

使用可能なリモートカーネルのリストです．

List of all available remote kernels

$Slaves

式の送受信
Sending and Receiving Expressions

基本的なリモート評価（ブロッキング）です．

Basic remote evaluation (blocking):

SetAttributes[RemoteEvaluate, HoldRest]

RemoteEvaluate[expr_, link_LinkObject]:=
    Module[{},
        Send[link, expr];
        Receive[link]
    ]

RemoteEvaluate[cmd]はこれをすべてのカーネルに送り，すべての結果を与えます．

RemoteEvaluate[cmd] sends it to all kernels, gives list of all results.

ここでカーネルをいくつか起動してみます．

Start a few kernels now...

TableForm[RemoteEvaluate[{$MachineName, $SystemID, $ProcessID, $Version}],TableHeadings -> {
            Automatic, {"host", "system", "process", 
              "Mathematica Version"}}]

プロセスのキュー（待ち行列）と並行性
Process Queues and Concurrency

ジョブをキューに追加，（いくつかが）終了するまでの待機，結果の取得

Queue jobs, wait for (some of) them to finish, get results.

ジョブ（プロセス）の開始と結果待ち
Starting Jobs (Processes), Waiting for Results

Queue[cmd]

与えられた式をキューに入れ，プロセスID(cookie)を返します．

Enters the given expression into queue, returns a process id (cookie).

WaitOne[pid_list]

与えられたプロセスの1つの終了を待ち， {result,pid,remaining_ids}を返します．

Waits for one of the given processes to finish and returns {result,pid,remaining_ids}.

Wait[pid_list], Wait[pid]

すべて，あるいは1つのプロセスを待ち，それぞれその結果（のリスト）を返します．

Wait for all processes or a single process, return list of results and single result, respectively.

連想配列は，様々な状態におけるすべてのプロセスの状態を監視します．まず，プロセスがキューに入り，それから特定のリンクで実行されます．さらに，終了し，結果の回収を待ちます．

Anassociative array keeps track of the state of all processes in the various states; first, a process is queued, then it is running on a certain link, then it is finished and waits for collection of result.

queued, link, finished[result]

ジョブをキューに入れます．

Enter a job into queue

SetAttributes[Queue, HoldAll]

Queue[cmd_] :=
  Module[{pid},
    AbortProtect[
      AppendTo[$Queue, Hold[cmd, pid]];
      qstate[pid] = queued ];
      trace[1, pid, " queued (", Length[$Queue], ")"];
    QueueRun[]
    pid
  ]

キューの管理とジョブのディスパッチ
Managing the Queue, Dispatching Jobs

リンクが使用可能になると，キューのプロセスを開始します．
結果が待機していると，それを保管し，使用中でないリンクのプールへリンクを戻します．

When link becomes available, start a queued process.
When result is waiting, store it and put link back on pool of idle links.

例
Examples

リモートホストが使用不能であれば，ローカルホストに2つのスレーブを作ります．

Create two slaves on local host if remote ones are not available.

slaves = 
      Table[LaunchSlave["localhost", "math -mathlink"], {2}]

ジョブをキューに入れ，終了を待ちます．

Queue a job and wait for it to finish.

j1 = Queue[Unevaluated[1 + 1]]

Wait[j1]

いくつかのジョブをキューに入れます．

Queue several jobs

Queue/@Hold[1 + 1, 2 + 2, 3 + 3]

これらを解放します．

Let them loose.

pids = List@@%

非決定的なジョブ待ちと，キューを徐々に空にするための戻り値の便利な使用法です．

Nondeterministic waiting for a job and a convenient use of the return value to drain queue:

{res, pid, pids} = WaitOne[pids]

クリーンアップします．

Cleanup:

Wait[pids]

レイテンシ隠蔽
Latency Hiding

レイテンシ：通信のオーバーヘッド，要求の処理の終了点と新しい要求の開始点の間の待ち時間．

Latency: communication overhead, the wait between end of servicing a request and start of servicing new one.

MathLink はバッファ付きのストリームです．これはすべての未処理の結果を受け取る前に，追加の要求を送ることができます．常にバッファを完全に満たしておきます．そのため遅延がありません．

MathLink is a buffered stream. Can send additional requests before receiving all outstanding results.
Keep buffers full at all times: no latency.

明示的なスケジューリング
Explicit Scheduling

1つのプロセッサにつき，1つのジョブを与えます．

One job per processor...

RemoteEvaluate[driver[]]

各プロセッサにプロセッサ番号を与えます．

Tell them their processor #.

ParallelMap[Queue[driver[#]]&, Range[Length[$Slaves]]]

並列評価
Parallel Evaluation

式の要素の並列評価
Evaluating Elements of an Expression in Parallel

ParallelEvaluate[h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]]

(Subscript[e, i])を使用可能なリンクに送り，結果を回収します．

Sends the (Subscript[e, i]) to available links, then collects results.

前の例の基本的な考え方を発展させました．

We developed the basic idea in a previous example.

SetAttributes[ParallelEvaluate, HoldAll]

ParallelEvaluate[h_[exprs___]] :=
      h @@ Wait[ List @@ Queue /@ Hold[exprs] ]

並列マッピング
Parallel Mapping

ParallelMap[f, h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]]

Map[f, h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]]を並列で評価します．

Evaluates Map[f, h[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]] in parallel.

これをParallelEvaluateに変えることもできますが，直接表現するのは非常に簡単です．

We could reduce it to ParallelEvaluate, but it's so simple to express directly:

ParallelMap[f_, h_[exprs___]] :=
      h @@ Wait[ Composition[Queue,f] /@ {exprs} ]

並列表
Parallel Tables

反復子変数の動的有効範囲を維持して，最も外側の次元を並列化します．

Parallelizing the outermost dimension, preserving dynamic scoping of the iterator variable

ParallelTable[e_, {i_Symbol, range__}, rest___List] :=
    ParallelMap[ Block[{i=#}, Table[e, rest]]&, Range[range] ]

例とトレース
Example and Traces

ParallelMap[Function[i, (i^2)], {1, 2, 3, 4}]

$TraceLevel = 1;

$TraceLevel = 0;

リモートカーネルの初期化
Initializing Remote Kernels

グローバル関数と変数の定義をスレーブに移動します．

Move definitions of global functions and variables to the slaves:

ExportEnvironment[{(Subscript[f, 1]), (Subscript[f, 2]), …, (Subscript[f, n])}]

ExportEnvironment["context`"]

Mathematica コードはポータブルです．スレーブにマシン固有コードをもつ必要がありません．
コードを転送するのに，既存のMathLink 接続を使います．分散の問題や共通FSの必要はありません．ツールキット自体のクライアントコードも，このリンクを介して転送されます．

Mathematica code is portable: no need to have machine-specific code on slaves!
Use existing MathLink connection to transport code: no distribution problems, no common FS needed. Even the client code for the toolkit itself is transported over this link.

例：行列の固有値
Example: Eigenvalues of Matrices

パラメータと2つの関数
A parameter and two functions.

prec = 18;

mat[n_] := Table[Random[Real, {-1, 1}, prec], {n}, {n}]

tf[n_] := Timing[Eigenvalues[mat[n]]][[1]]

すべてのスレーブに知らせる
Make them known to all slaves.

ExportEnvironment[{prec, tf, mat}];

確認します．

Check:

Clear[prec];

RemoteEvaluate[prec]

試行
A sample run

ParallelMap[tf, {5, 10, 15, 20, 25}]

例：無限大の計算
Example: an Infinite Calculation

これらの多項式は約分できないと思われます．

These polynomials are conjectured to be irreducible.

poly[n_, x_] := (∑_(i = 1)^(n + 1))i (x^(i - 1))

ExportEnvironment[{poly}];

$LoadFactor = 2;

明示的にスケジューリングします．すべてのプロセッサが使用中になるようにします．しかし，長い入力キューは避けます．

Explicit scheduling; try to keep all processors busy, but avoid long input queues.

d = 150; ids = {};
While[True,
\(While[$QueueLength == 0, 
	
                With[{d = d}, 
                  AppendTo[ids, 
                    Queue[{d, 
                        Length[
                          FactorList[
                            poly[d, 
                              x]]]}]]
]; 
\(d++\)
]\);
If[
                Length[ids] > 0, 
{res, id, ids} = 
                  WaitOne[ids]; 
Print[
                  res]
];
	]

割込みの後はオーファンを回収します．

Collect orphans after an interrupt:

Wait[ids]

仮想共有メモリ
Virtual Shared Memory

プロセスは互いにどのように通信するのでしょうか．

How can processes communicate with each other?

仮想共有メモリ：分散メモリマシンに擬似的に作成．

Virtual Shared Memory: faking it on distributed memory machines.

プロセスがグローバル（共有）変数の値を必要としたり，値を代入したりする場合，プロセスはスケジューラに要求を送り返します．スケジューラは変数を管理し，値を更新し，要求された値を送り返します．

When a process needs the value of a global (shared) variable or wants to set its value, it sends a request back to the scheduler. The scheduler maintains the variables, udpates their values, and sends back the requested values.

共有変数・関数の宣言
Declaring Shared Variables/Functions

宣言により，スレーブ側の定義が設定され，すべてのスレーブに転送されます．従って，別の初期化コードは必要ありません．

A declaration sets up the slave-side definitions and transports them to all slaves, which, therefore, need no separate init code.

SharedVariables[{(Subscript[v, 1]), (Subscript[v, 2]), …, (Subscript[v, n])}]

SharedDownValues[{(Subscript[f, 1]), (Subscript[f, 2]), …, (Subscript[f, n])}]

例
Examples

ここでカーネルを再び初期化します．

Re-init kernel here...

Needs["Parallel`VirtualShared`"]

共有変数を宣言します．

Declare shared variables.

SharedVariables[{x, y, z}]

基本的なテストを行います．

Basic tests.

x = 0;

Eval[links[[1]], x]

Eval[links[[1]], x = x + 1]

x

原始的な更新

Atomic update:

Eval[links[[2]], x++]

x

共有変数を使うのに，特別な文法は必要ありません．通常の方法で行います．変数が共有されると，すべてのアクセスがマスタープロセスに透過的に送られます．

No special syntax is required to use shared variables. Just do things the normal way. If a variable is shared, all accesses are transparently sent to the master process.

原型的同期問題：読取りと書込みの間に，その他のプロセスが値を変更したり，古い値を読んだりすることが可能．

The prototypical synchronization problem: between reading and writing, another process can change the value or read the obsolete value.

y = 0;
ParallelMap[(
		Pause[0.1 Random[]];
		a=y;
		Pause[Random[]];
		y=a + 1)&,
	Range[10]]

y

同期
Synchronization

変数への排他的なアクセスは新しい基本処理，「conditional atomic write」で行われます．機械語では，「test and set」あるいは「read-modify-write」と呼ばれることもあります．

Exclusive access to variables can be achieved with one new primitive operation, a conditional atomic write. In machine code, it is often called test and set or read-modify-write.

TestAndSet[s, e]

現在sの設定が何もない場合，あるいは値がNullの場合のみ，（s=eのように）sにeを代入します．sの（新しい，あるいは古い）値を返します．read-modify-writeは原始的です．

Set s to e (just like s=e), but only if s is currently unset or its value is Null. Return the (old or new) value of s. The read-modify-write is atomic.

l1 = $Slaves[[1]];  l2 = $Slaves[[2]];

SharedVariables[{z}]

Eval[l1, TestAndSet[z, 17]]

Eval[l2, TestAndSet[z, 18]]

Eval[l1, z = Null]

Eval[l2, TestAndSet[z, 18]]

重要なセクション
Critical Sections

変数やその他のリソースへの排他的アクセスは「lock」で制御できます．以下に，変数への非同期アクセスについての前の例に追加される典型的なコードを示します．

Exclusive access to a variable or other resource can be controlled by a lock. Here is the typical code added to the previous example about asynchronous access to a variable.

SharedVariables[{y, lock}]

$TraceLevel = 2;

y = 0;
ParallelMap[(
		Pause[0.2  Random[]];
		While[TestAndSet[lock, #] =!= #, Pause[0.8 Random[]]];
		a = y;
		Pause[Random[]];
		y = a + 1;
		lock = Null;
		a + 1)&,
	Range[5]]

y

$TraceLevel = 0;

食事をする哲学者
The Dining Philosophers

共有リソースへのアクセスの古典的な例として，テーブルの周りに座っている哲学者たちの間に置かれたフォークが挙げられます．1人の哲学者につき，左手と右手で食べるために2本のフォークが必要です．

The classic example for access to shared resources, in this case forks placed between philosophers sitting at a round table. A philosopher needs the two forks to his left and right to eat.

変数forks[1],…,forks[n]は，それらを条件付きで固有の値に設定することで，フォークを取っておく役割を果たします．

The variables forks[1],…,forks[n] serve to reserve forks by conditionally setting them to a unique value.

SharedDownValues[{forks}]

stop = False; SharedVariables[{stop}]

それぞれの哲学者が1つのプロセスを示します．

Each dining philosopher is one process.

n = Length[$Slaves]

以下に哲学者の考え方と食べ方に関する動作を示します．

Here's the behaviour of a philosopher: either thinking or eating.

philo[i_]:=
With[{f1 = i, f2 = Mod[i, n] + 1},
While[! stop,
Print[i, " thinking"];
  
Pause[10 Random[]];
	 	  
Print[i, " hungry"];
			While[TestAndSet[forks[f1], i] =!= i || 
                        TestAndSet[forks[f2], i] =!= i,
Print[i, " angry"];
Pause[Random[]]];
  
Print[i, " eating"];
Pause[2 Random[]];
  
forks[f1] = Null;forks[f2] = Null;
]]

ExportEnvironment[{philo, n}];

本質的にはParallelMapが必要ですが，これらを手際よく止める方法が必要です．

We want essentially ParallelMap, but need a way to stop them cleanly.

ParallelMap[philo, Range[n]]

pids = Composition[Queue, philo] /@ Range[n]; Wait[pids]

stop = True; Wait[pids]

並列プログラミング
Parallel Programming

ここでは，並列Forループ以上のことを取り上げます．

More than parallel For loops.

プロセス生成の関数スタイル
A Functional Style for Process Generation

例えば，一般化された内積が必要であるとします．

Say, you want a generalized inner product,

(Subscript[((a . b)), (Subscript[i, 1]) (Subscript[i, 2]) …(Subscript[i, n - 1]) (Subscript[k, 2]) 
          (Subscript[k, 3]) …(Subscript[k, m])]) = 
      p[t[(Subscript[a, (Subscript[i, 1]) (Subscript[i, 2]) …(Subscript[i, n - 1]) 1]), 
          (Subscript[b, 1  (Subscript[k, 2]) (Subscript[k, 3]) …(Subscript[k, m])])], …, 
        t[(Subscript[a, (Subscript[i, 1]) (Subscript[i, 2]) …(Subscript[i, n - 1]) d]), 
          (Subscript[b, (Subscript[dk, 2]) (Subscript[k, 3]) …(Subscript[k, m])])]]

ここで，dはaの共通の最終次元であり，bの最初の次元であるとします（pを和，tを積とお考えください）．

where d is the common last dimension of a and first dimension of b.
(Think of p as plus and t as times.)

Inner[t, Array[a, {2, 2}], Array[b, {2, 2}],  p]

並行実行のためp[…]のすべての計算をキューに入れるために，pの代りにComposition[Queue,p]を用いることもできます．その結果，プロセスIDのテンソルが得られます．

You can use Composition[Queue,p] in place of p to cause all calculcations of p[…] to be queued for concurrent execution.
The result is a tensor of process IDs.

skel = 
      Inner[t, Array[a, {2, 2}], Array[b, {2, 2}], Composition[Queue, p]]

プロセスIDをリストに集め，すべてのプロセスが関連の結果のリストを返すのを待ちます．

Collect the pids in a list, and wait for all processes, returning the list of corresponding results.

pids  =  Flatten[skel]

ress  =  Wait[pids]

最後に，置換を用いてテンソルのプロセスIDを結果と置き換えます．

Finally, replace the pids in the tensor by the results, using substitution.

skel  /.  Thread[pids -> ress]

このスケジューリングは，通信に関しては最適ではありません．よりよい可能性としては，aの次元に沿ってのみ並列化し，テンソルbをジョブのスケジューリングの前に各プロセッサに送るということが挙げられます．

This scheduling is not optimal in terms of communication. One better possibility is to parallelize only along the dimensions of a and send the tensor b once to each processor before scheduling the jobs.

明示的なスケジューリング：Gaxpy
Explicit Scheduling: Gaxpy

「Golub&van Loan, Matrix Computations, (3^rd)ed., The Johns Hopkins Univ. Press, 1996」より抜粋．

Taken verbatim from Golub&van Loan, Matrix Computations, (3^rd)ed., The Johns Hopkins Univ. Press, 1996.

Y=Y+A.X, A an n×n matrix, X, Y, n-vectors

Needs["Parallel`VirtualShared`"]

サンプル（記号）行列とベクトルです．

Sample (symbolic) matrices and vectors.

n = 4;

A = Array[a, {n, n}];
X = Array[x, {n}];
Y = Array[y, {n}];

SharedVariables[A, X, Y];

スケジューリングします．nの約数であるプロセッサの数を求めます．

Scheduling. Find a number of processors that is a divisor of n.

\(p0 = Length[$Slaves]; 

    p = Max[Select[Divisors[n], # ≤ p0&]]\)

r = n/p

クライアントコードです．

The client code

\(SetAttributes[gaxpy, HoldAll]; 

    gaxpy[μ_, A_, X_, Y_] := 
	
      Module[{Xl, Yl, j, Al}, 
		
        With[{row = Range[\((μ - 1)\) r + 1, μ\ r]}, 
		Xl = X; 

          		Yl = Y[[row]]; 
		
          Do[
			Al = A[[row, j]]; 

            			\(Yl = 
              Yl + Al*Xl[[j]]; \)

            			, {j, 1, n}]; 
		
          Y[[row]] = Yl; 
		]]\)

コード（データではありません）をエキスポートします．

Export the code (not the data!).

ExportEnvironment[p, n, r, gaxpy];

すべてのリモートプロセスを開始します．

Start all remote processes.

Wait[Queue[gaxpy[#, A, X, Y]] & /@ Range[p]];
Y

$TraceLevel = 1;

$TraceLevel = 0;

処理能力測定
Performance Measurements

オーバーヘッド
Overhead

普通の1プロセスの処理時間を求めます．

Time for servicing one trivial process.

With[{n = 200, t = AbsoluteTime[]}, 
ParallelMap[
        Identity, 
        Range[n]]; 
((AbsoluteTime[] - t)/n)]

スピードアップと効率
Speedup and Efficiency

スピードアップとは，全並列実行時間の順次実行時間に対する比です．

Speedup is the ratio of total parallel execution time to sequential time.

効率とは，実際のスピードアップのプロセッサ数に対する比です．

Efficiency is the ratio of observed speedup to number of processors.

このコードは普通の個々の並列計算のスピードアップと効率を測定します．

This code measures the speedup and efficiency for a trivial independent parallel computation.

Performance[ntry_, d_:  1]:=
Module[{p = Length[$Slaves], t, su},
t = AbsoluteTime[];
ParallelMap[Pause[d] &,Range[ntry]];
su = ((ntry d)/(AbsoluteTime[] - t));
{p, su, su/p}]

Performance[12, 4]

アプリケーション例
Application Examples

最小表面

A minimal surface

並列アニメーション

A parallel animation

おわりに
Conclusions

これで，Mathematica 内部だけでインタラクティブな並列記号，数値，画像計算ができるようになり，大変うれしく思います．

I am excited that one can now do interactive parallel symbolic, numeric and graphic computation entirely within Mathematica.

このツールキットにより，広範囲のワークグループ，研究室，教室での本格的な並列計算の利用が簡単になります．これは，基本的には既存のMathematica ユーザがすぐに使えるものです．

This toolkit should make serious parallel computing accessible to a wide range of workgroups, labs and classrooms." It is basically ready for all existing Mathematica users.