LU分解アルゴリズムのおさらい

LU分解を使って連立一次方程式を解く方法

LU分解は係数行列Aを、A=LUと、上三角行列Lと下三角行列Uに分解して解くものでした。

例えば係数行列Aのサイズが3×3の場合は

LUeq02

と分解できます。

これを使って連立一次方程式

……式02

を解くには、

……式03

となるベクトルcをおいて、式02に代入し、

……式04

より、ベクトルcの値を求めます。この処理を前進代入といいます。

得られたベクトルcの値を式03に代入すると変数ベクトルxの値が求まります。この処理を後進代入といいます。

LU分解の方法

LU分解の方法はいくつかあります。今回は片桐さんの資料でも並列化に向くと紹介されている、外積形式ガウス法を用いて計算を行います。

外積形式ガウス法は英語ではDoolittle algorithmと呼ばれています。そのアルゴリズムは、ガウスの消去法と同様の操作でLU分解を行うというものです。以下でそのアルゴリズムを説明します。

n行n列の係数行列A、変数ベクトルx、定数項bを与えられ、これを解いていきます。

ガウスの消去法の前進消去は（n-1）個の変数を、1度の処理で1つずつAから削除していきます。この削除の処理を何度行ったかわかりやすくするため、削除処理を行った回数をAの上につけて表します。初期状態では次のように定義されます。

……式05

次にA^(k-1)のとき、k行ｋ列の変数を削除する処理を考えます。

eq07 ……式06

この削除処理は対角線上の係数a_kkを用いて列方向に行います。この消去を列方向のインデックスi(i=k+1, k+2, …, n)を用いて表すと式07になります。

……式07

ここからは式07の処理を、行列で表していきます。

行列Lを式08と定義すると、A^(k)に対する消去は式09のように表せます。

eq09 ……式08

……式09

この式09より、式07を行列式で表せました。

次に係数行列Aをガウスの消去法で削除していく全過程を式09を使って表します。つまり初期状態の A^(0) から A^(n-1) までの消去の式を式09を用いて表すと、次の式10となります。

eq11 ……式10

この式10が終了する――つまりガウスの消去法の前進消去が終了すると、上三角行列Uが得られます。つまりは、

……式11

となり、上三角行列Uは求まります。

下三角行列Lは、逆行列L^(-1)の要素の符号を変転させたものなので、容易にもとまります。

以上で、LU分解のアルゴリズムの説明を終わります。

MPIを用いた並列プログラミングの概要

ここでは並列計算機の種類とデータ通信の方法から、 MPIを用いた並列プログラミングの概要について説明します。

並列計算機の種類とMPI

並列計算機には共有メモリ型と分散メモリ型の2つの型があります。前者の共有メモリ型は、すべてのプロセッサが同じメモリを参照する構成の計算機です。後者の分散メモリ型は、各CPUがそれぞれメモリを持っており、各CPUがお互いのメモリを直接参照できない構成の計算機です。

今回製作したRaspberry Piスパコンは分散メモリ型の並列計算機です。分散メモリ型の並列計算機が各CPUのメモリの中身を参照するには、メモリの中身をメッセージという形で交換しあう必要があります。このメッセージを交換することをメッセージ・パッシングと呼び、これを行う規格の１つがMPIです。

MPIはあくまで規格であり、その実装は多々あります。その実装されたライブラリの一つが、PHASE0やLinpackで使用していたMPICHです。

MPICHの最新はMPICH3ですが、 Raspbianに提供されているバージョンはMPICH2（v1.4.1）なので、今回はMPICH2を用います。

MPI関数の概要

MPIには数百種の関数がありますが、そのすべてを覚えるのは大変です。しかし、大抵の並列プログラムは数種類の関数を知っていれば実装できるので、ここからはその最低限の関数について説明します。

すべての関数の詳細が知りたい方は、MPICHのドキュメントを参照してください。

http://www.mpich.org/static/docs/v3.1.3/

システム関数

MPIを利用するための関数です。初期化を行うMPI_Init関数や、並列処理を行うプロセス全体の数を調べるMPI_Comm_size関数、プロセス全体のうち、自分のプロセスに付けられた番号を調べるMPI_Comm_rank関数、終了処理を行うMPI_Finalizeなどがあります。

1対1通信関数

MPIでの通信はプロセス単位で行われます。プロセスとはOSよりメモリが与えられて実行されているプログラムのことです。今回使用するRaspberry Piクラスタでは、1ノードにつき1プロセス立ち上がっています。そのため、プロセス間の通信を行う説明をした場合は、異なるRaspberry Pi間での通信が行われていると考えてください。

このプロセス間を1対1で通信する関数を1対1通信関数と呼びます。大別すると以下の2種があります。

ブロッキング型関数

データの送受信が完了するまで、1対1通信関数を呼び出した箇所で処理を停止する関数です。送信はMPI_Send、受信はMPI_Recvなどが該当します。

ノンブロッキング型関数

データの送受信が完了するのを待たずに、すぐに1対1通信関数を呼び出した箇所に戻る関数です。データの送受信が完了したかはMPI_Wait関数を用いて確認します。送信はMPI_Isend、受信はMPI_Irecvなどが該当します。

1対全通信関数

1つのプロセスからのメッセージを、並列処理全体を構成する全プロセスへ放送する関数です。 1対全通信関数は放送関数とも呼ばれます。 MPI_Bcast関数が該当します。

MPIを用いた並列HelloWorld

具体的なMPI関数の使い方を説明するため MPI版並列HelloWorldを作りました。以下にプログラムを示します。

#include 
#include "mpi.h"

int  main(int argc, char* argv[]) {
    int myid, numprocs;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

    printf("Hello world! Myid: %2d / %2d \n", myid, numprocs);

    MPI_Finalize();
}

ソースはこちらからもダウンロードできます。
mpi_hello.c

MPIの初期化と終了

MPI関数は、初期化を行うMPI_Init関数と、終了処理を行うMPI_Finalize関数の間でのみ使用できます。

int main(int argc, char* argv[]){
    MPI_Init(&argc, &argv);

    // ここにコードを書く

    MPI_Finalize();
}

MPI_Init関数の引数には、main関数の引数argcとargvのアドレスを入れます。これによってMPI_Init関数は、起動時のコマンドラインオプションを読み込み、いくつかの設定を行うことができます。しかし今回は特に必要ないので、気にしないで大丈夫です。

MPI_Finalize関数はMPIの終了処理を行い、MPIを停止します。

プロセス番号の取得

mpirunコマンドでMPIのプログラムが実行されると、並列計算を行うすべてのノード上で「同じプログラム」が一斉に起動します。例えば、MPI版並列HelloWorldを16並列で実行した場合、 16台すべてのマシンでMPI版並列HelloWorldが起動します。

起動した各プロセスにはそれぞれ「ランク」という値を割り当てられ、これを取得するMPI関数がMPI_Comm_rank関数です。ランク値は自身に割り当てられた仕事範囲の確認や、他のプロセスを指定しての通信を行うために使われます。

先ほどのHelloWorldでは、MPI_Comm_rank関数は次のように用いられています。

MPI_Comm_rank(MPI_COMM_WORLD, &myid);

第1引数にはコミュニケーターという通信集団を指定します。コミュニケーターを分けることで高度な通信が可能となりますが、小規模な計算では必要ありません。 MPI_COMM_WORLDを設定して、すべてのプロセスの所属する通信集団を指定します。

第2引数にはランク値を保存する変数のアドレスを指定します。今回の場合、変数myidに自プロセスのランク値が設定されます。

総プロセス数の取得

MPI_Comm_size関数はコミュニケーターに所属するプロセスの数を取得します。先ほどのHelloWorldでは、次のように用いられています。

MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

第1引数はMPI_Comm_rank関数と同様にコミュニケーターの指定を行います。

第2引数にはプロセス数を保存する変数のアドレスを指定します。今回の場合、変数numprocsに全プロセス数が設定されます。

printfの出力先

mpirunで起動したプロセス内でprintfによる出力を行うと、その出力はプロセスを実行している各ノード上ではなく、 mpirunを実行したノードの標準出力へ出力されます。

よってMPI版HelloWorldの次のprintf文は、各プロセス上で取得された自身のランク値を出力し、その出力はmpirunを実行したノードの標準出力へ出力されます。

printf("Hello world! Myid: %2d / %2d \n", myid, numprocs);

コンパイルと実行方法

先ほどのMPI版並列HelloWorldコードをRaspberry Pi上の ~/mpi_hello/mpi_hello.c に保存し、ディレクトリを移動してください。

$ cd ~/mpi_hello/

MPICH2でMPIを用いたプログラムのコンパイルは、以下のコマンドで行います。

$ mpicc -o mpi_hello mpi_hello.c

この作業を並列計算を行うマシンで行った後、以下のコマンドを実行すると、MPI版HelloWorldが16ノードで動作します。

※第2回で作成したmpdhostsファイルが ~/mpdhosts に保存されているとします

$ mpirun -f ~/mpdhosts -np=16 ~/mpi_hello/mpi_hello

出力結果は以下のようになり、ランク値の順番が一部バラバラになっていることからもプログラムが並列に動作したことが感じられます。

Hello world! Myid:  0 / 16
Hello world! Myid:  1 / 16
Hello world! Myid:  2 / 16
Hello world! Myid:  3 / 16
Hello world! Myid:  5 / 16
Hello world! Myid:  4 / 16
Hello world! Myid:  6 / 16
Hello world! Myid:  7 / 16
Hello world! Myid:  9 / 16
Hello world! Myid:  8 / 16
Hello world! Myid: 11 / 16
Hello world! Myid: 12 / 16
Hello world! Myid: 10 / 16
Hello world! Myid: 14 / 16
Hello world! Myid: 13 / 16
Hello world! Myid: 15 / 16

MPI関数を用いたデータ通信サンプルプログラム

通信関数を用いたサンプルを作って、通信関数の使い方を学びます。

1対1通信関数の説明

MPI_SendとMPI_Recvは、プロセス間を1対1で通信する関数です。

最初にMPI_Sendの関数定義を示します。

int MPI_Send(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm)

引数は以下になります。

第1引数: 送信するデータのアドレス
第2引数: 送信するデータの個数
第3引数: データタイプ
第4引数: 送信先のランク値
第5引数: メッセージタグ（今回は説明を省略。0を設定する）
第6引数: コミュニケータ
戻り値: エラーコード

第3引数のデータタイプには、MPI_CHAR, MPI_INT, MPI_FLOAT, MPI_DOUBLEなどを指定します。これらは基本的に、C言語のchar, int, float, doubleに対応しています。

次にMPI_Recvの関数定義を示します。

int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status)

引数は以下になります。

第1引数: 受信先のアドレス
第2引数: 受信するデータの個数
第3引数: データタイプ
第4引数: 送信元のランク値
第5引数: メッセージタグ（今回は説明を省略。0を設定する）
第6引数: コミュニケータ
第7引数: 受信状況に関する情報
戻り値: エラーコード

基本はMPI_Sendと同様です。ひとつ違うのは第7引数で、ここで指定したMPI_Status型の変数に受信状況に関する情報が入ります。詳しくは説明を省略します。

放送関数の説明

MPI_Bcast関数は、特定のプロセスからその他すべてのプロセスに、データを送信する関数です。

以下にMPI_MPI_Bcastの関数定義を示します。

int MPI_Bcast( void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm )

引数は以下になります。

第1引数: 送信または受信するデータの先頭アドレス
第2引数: (送信または受信する)データの個数
第3引数: データのタイプ
第4引数: 送信元プロセスのランク値
第5引数: コミュニケーター

MPI_Bcast関数は、送信も受信も同じ関数で行います。送信受信の役割を分けているのは、第4引数の値です。ここで指定したランク値のプロセスから送信が行われ、他のプロセスは受信を行います。
第1引数には、送信する側（第4引数で指定するランク値のプロセス）は送信したいデータの開始アドレスを、受信する側は受信先の開始アドレスを設定します。

サンプルプログラム

以上に説明した通信関数を用いて、変数aの値をやりとりするサンプルを作りました。

#include 
#include "mpi.h"

int  main(int argc, char* argv[]) {
    int myid, numprocs;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

    int a = 10;

    if(myid == 0){
        // プロセス0だけaに0を代入
        a = 0;
    }

    // プロセス0からaの値を放送
    MPI_Bcast( &a, 1, MPI_INT, 0, MPI_COMM_WORLD );

    // ここで放送を受けた全プロセスのaの値は0になる
    printf("a = %d, Myid: %2d / %2d \n", a, myid, numprocs);


    // aにランク値を追加
    a += myid;

    // 1対1通信関数でプロセス1のaの値をプロセス0に送信する
    if(myid == 0){
        MPI_Status istatus;
        MPI_Recv(&a, 1, MPI_INT, 1, 0, MPI_COMM_WORLD, &istatus);
        // プロセス1から受け取ったaの値を表示
        printf("Receive from proc 1\n");
        printf("a = %d\n",a);
    }
    if(myid == 1){
        MPI_Send(&a, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
    }

    MPI_Finalize();
}

以下、コードの説明をします。

最初に変数aを定義します。このaの値は各プロセスにより異なり、プロセス0はa=0、それ以外はa=10となります。

次に放送関数を用いて、プロセス0のaの値をその他すべてのプロセスに放送します。放送を受けたプロセスは、aの値が書き換えられ、プロセス0と同じa=0となります。

次に変数aにプロセスのランク値を加算した後、1対1通信関数を用いて、変数aの値をプロセス1からプロセス0に送ります。これにより、プロセス0の変数aの値は0から1に書き換えられます。

このプログラムを3ノードで実行すると、次のような結果が得られます。

$ mpiexec -np 3 -f ../mpdhosts ./mpi_communication
a = 0, Myid:  2 /  3
a = 0, Myid:  0 /  3
a = 0, Myid:  1 /  3
Receive from proc 1
a = 1

以上でMPIとMPIを用いた並列プログラミングについての、簡単な解説を終わります。

RaspberryPiクラスタ製作記第3回「並列プログラミング」

第３回　並列計算プログラムをつくろう

前回までは既存のプログラムを使い、並列化による性能上昇について調べました。これにより大雑把に「並列化するとプログラムの処理を高速化できる」ことはわかりましたが、具体的にどのようにすれば自分のプログラムを並列化して高速化できるのかは、まだわかりません。

そこで今回からは、連立一次方程式を解くプログラムの制作を通して、並列計算プログラムを1から自作する方法や、高速化のための手法について学んでいきたいと思います。

なぜ連立一次方程式を扱うかというと、それは理解が容易だからです。スパコンでは様々な計算が行われていますが、大抵これらの計算は非常に難しく、短時間で理解し、プログラムを作るのは大変です。しかし、連立一次方程式は小さいサイズなら筆算で簡単に解けるため、これを解くプログラムの理解も容易です。加えて、連立一次方程式をプログラムで解く方法は、数値計算の教科書や並列計算の入門本など、非常に多くの書籍で扱われています。そのため、内容でわからないことが出た場合や、より詳しく知りたい場合などは、豊富な参考文献をあたって調べられます。

本記事も、片桐孝洋さんの『スパコンプログラミング入門並列処理とMPIの学習』という本と、同じく片桐さんの講義資料を参考にしています。わからないことが出た際はこちらを参照してください。

LU分解法（１）
LU分解法（２）
LU分解法（３）
東京大学情報基盤センター准教授片桐孝洋

連立一次方程式をプログラムで解く

今回は「LU分解」を用いて、以下の簡単な連立一次方程式を解くプログラムを作っていきます

一般に、問題を解くプログラムのことをソルバと呼び、今回の連立一次方程式を解くプログラムのことは連立一次方程式ソルバと呼びます。

LU分解のアルゴリズムはいくつかあります。今回は片桐さんの資料でも並列化に向くと紹介されている、外積形式ガウス法を用いて計算を行います。

※外積形式ガウス法を用いたLU分解のアルゴリズムについては下記のページを参照してください
アルゴリズムのおさらい

行列の定義

係数行列A、変数ベクトルx、定数項b、ベクトルcはプログラムで次のように定義します。

// SIZE: 係数行列Aのサイズ
#define SIZE 3

double A[SIZE][SIZE] = {
    {4, 1, 1},
    {1, 3, 1},
    {2, 1, 5}
};

double b[SIZE]={9,10,19};

double c[SIZE];
double x[SIZE];

LU分解のコードを作る

最初に係数行列Aを下三角行列Lと上三角行列Uに分解するコードを作ります。

下三角行列Lと上三角行列Uは一般的に1つの行列で表されます。

行列LUの要素は、LU分解より係数行列Aを更新しつつ作られるので、係数行列Aに上書きして入れていきます。そして、LU分解が完了すると、行列Aの中身はすべて行列LUの要素におきかわります。その後の前進代入と後進代入は、行列Aを置き換えて作られた行列LUを用いて行います。

外積形式ガウス法を用いたLU分解は、行方向のループk（k=0,1,…,n=SIZE)と、そこに入れ子になった以下の2つのループで行います。

・行列L（枢軸ベクトル）を作っていく列方向のループ……(1)
・行列Uを作り、残りの値を更新する行と列の2重ループ……(2)

です。

この時、1つ目のループで計算されるベクトルのことを枢軸(すうじく)ベクトルと呼びます。
以上の処理をプログラムで表すと、次のようになります。

int n = SIZE;

// LU分解開始
for (k=0; k<n; k++) {
    // kで回す行方向ループ

    // (1)
    // 行列Lの処理
    double div_tmp = 1.0 / A[k][k];
    for (i=k+1; i<n; i++) {
        // iで回す行方向ループ

        // l_ik = a_ik / a_kk
        A[i][k] = A[i][k] * div_tmp;
    }

    // (2)
    // 行列Uの処理
    for (j=k+1; j<n; j++) {
        // jで回す行方向のループ

        // Tips: ここで値を変数に入れたほうが
        // ループ内で参照するよりメモリアクセス回数が減ってお得
        double a_jk = dtemp = A[j][k];

        for (i=k+1; i<n; i++) {
            // iで回す列方向のループ
            A[j][i] = A[j][i] - A[k][i] * a_jk;
        }
    }
}
// LU分解終了

前進代入

LU分解の後は下三角行列Lと定数項bの値より、Lc=bのベクトルcの値を求めます。

// 前進代入
// Lc = b
for(k=0; k<n; k++){
    // 行方向のループ
    c[k] = b[k];
    for(i=0; i<k; i++){
        // 列方向のループ
        c[k] -= A[k][i] * c[i];
    }
}

後進代入

前進代入で求めたベクトルcより、Ux = cを解いて変数ベクトルxを求めます。

// 後進代入
// Ux = c
// k = n - 1
x[n-1] = c[n-1] / A[n-1][n-1];
// k = (n - 2) ~ (0)
for(k = n - 2; k >= 0; k--){
    x[k] = c[k];
    for(i=k+1; i<n; i++){
        // 列方向のループ
        x[k] -= A[k][i] * x[i];
    }
    x[k] = x[k] / A[k][k];
}

完成

以上でLU分解を用いた連立一次方程式ソルバは完成です。以下にコードを示します。

LU_normal.c

このコードを実行すると、式eq01の係数行列AをLU分解した結果と、変数ベクトルxの解が出力されます。

$ gcc -o LU_normal LU_normal.c

$ ./LU_normal
LU:
 4.00  1.00  1.00
 0.25  2.75  0.75
 0.50  0.18  4.36

x:
 1.00  2.00  3.00

LU分解コードを並列化しよう

ここからはMPIを用いて、連立一次方程式を並列に解くプログラムを作っていきます。ただし全ての並列化を行うのは大変なので、計算に一番時間がかかるLU分解の部分のみを並列化します。

※MPIについてご存じない方は、以下のページを参照してください。
MPIを用いた並列プログラミングの概要

最初に完成したコードを示し、このコードを扱いながら説明を行います。以下からダウンロードして下さい。

LU_mpi_3x3

担当範囲の計算

外積形式ガウス法を用いたLU分解の並列化を、今回は列方向に分割して行います。

各プロセスはそれぞれ、以下のプログラムで計算した、開始（start）から終了（end）までの列を担当し、計算します。担当列の長さは（partition_size）で表します。

// 1つのプロセスが担当する列の長さ
int partition_size = n / numprocs;

// プロセスが担当する列の開始点
int start = myid * partition_size;

// プロセスが担当する列の終了点
int end = (myid + 1) * partition_size;

行列Lの計算と放送

下三角行列L（枢軸ベクトル）の計算は、その列を担当するプロセスが行います。

枢軸ベクトルは、その他のプロセスが係数行列Aの更新に必要とするので、別のプロセスに放送します。そのため、枢軸ベクトルを担当したプロセスは、別途放送用の配列Lを用意し、更新した値を配列Lに入れて放送します。

この処理を行うのが、次のコードです。

int worker = k / partition_size;

// 行列Lの処理
if(worker == myid){
    // kが自分の担当範囲内の場合はLを処理
    double div_tmp = 1.0 / A[k][k];
    for (i=k+1; i<n; i++) {
        // iで回す行方向ループ

        // l_ik = a_ik / a_kk
        A[i][k] = A[i][k] * div_tmp;
        L[i] = A[i][k];
    }
}

// Lを放送
MPI_Bcast(&L[k], n - k, MPI_DOUBLE, worker, MPI_COMM_WORLD);

行列Uを計算する

枢軸ベクトルを受け取った各プロセスは、行列Uの計算と、行列Aの値の更新を行います。

最初に行うのは、計算に使用する枢軸ベクトルの要素a_jkの取り出しです。

double a_jk = L[j];

次に自身の担当する列のみ計算するため、i-ループを次のように書き換えます。

if(k < start){
    i = start;
}
else{
    i = k+1;
}

for (; i < end; i++) {
// 行列Uの計算と更新処理

kがプロセスの担当列番号より小さい場合は、startからendまでループを回して計算します。

その他の場合はkが自分の担当範囲か、自分の担当範囲が既に終了しているので、k+1からendまでループを回して計算します。

以上でLU分解の処理は終了です。

結果を集計する

今回の並列化はLU分解のみとしているので、 LU分解後の処理はすべてプロセス0が行います。そのため、LU分解後は各プロセスの計算した配列LUのデータを、プロセス0に集めます。

これを行うのが次のコードです。

/ 結果を集める
for(i = 0; i < n; i++){
    // 行方向のループ
    for(j = partition_size; j < n; j+=partition_size){
        // 列方向のループ

        int worker = j / partition_size;

        if(myid == 0){
            // process0 はその他のプロセスの担当部分を集める
            MPI_Recv(&A[i][j], partition_size, MPI_DOUBLE, worker, 0, MPI_COMM_WORLD, &istatus);
        }
        if(myid == worker){
            // process0 以外は、自分が担当した部分を送る
            MPI_Send(&A[i][j], partition_size, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
        }
        // すべてのプロセスがここに来るのを待つ
        MPI_Barrier(MPI_COMM_WORLD);
    }
}

性能測定

プログラムの並列化は少々大変でしたが、このプログラムを複数ノードで並列に動かすことで、高速に計算できるはずです。

係数行列のサイズは3程度では変化がわからないので、サイズを192に拡大して実行しました。その結果をグラフに示します。

PHASE/0のときのようにノード数が増えるごとに計算時間が減って……いませんね。今回のプログラムはPHASE/0の時と真逆で、ノード数を増やすごとに処理時間が増えています。

どうやらLU分解を並列化しただけでは、並列化のコスト（主に通信のコスト）のほうが大きいため、逆に時間が増えてしまうようです。

残りの前進代入と後進代入も並列化すれば、ノード数に応じて処理時間が減っていく、期待通りの結果が得られるかもしれません。余力のある方は、最初に挙げた片桐さんの資料を参考に、並列化を行ってみてください。私は別の方法で高速化を図ってみます。

BareMetalで作る連立一次方程式ソルバ専用機

現代のOSは様々な機能を私たちに提供してくれています。マルチタスクシステムもその１つです。このシステムは同じマシンを複数の人が使うことを想定し、プロセスを人間が分からない程度の短時間で切り替え、あたかも複数のプログラムが同時に動いているかのように実行するシステムです。これのおかげで私たちは先程のLU分解のプログラムを動かしている間も、別マシンからSSHアクセスして作業することが可能になっています。

しかし、今回のように私1人が1つのプログラムだけを処理したい場合、マルチタスクシステムを含め、OSの便利な機能は必要のない、無駄な機能です。このマルチタスク処理を含め、OSの実行にかかる様々なコストを削減して高速化できないでしょうか。

私はこの目的を達成する方法に、心当たりがあります。著者紹介欄でもご紹介頂いているように、私は「Baremetalで遊ぶ Raspberry Pi」という電子書籍を書きました。

通常Raspberry Piは、電源投入後にSDカード上のプログラムを読み込み、 LinuxなどのOS(カーネル)を起動します。ここで、SDカード上の設定を少し書き換えると、本来カーネルが起動するところを、代わりに自作のコードを実行できます。そうすればOSの実行にかかる様々なコストを削減して、計算を高速に行えるはずです。

この自作のコードの作り方と、OSを用いない素のハードウェア（BareMetal）上で、 Raspberry Piを制御する方法を説明した電子書籍が「Baremetalで遊ぶ Raspberry Pi」です。

ここからは前半に作った連立一次方程式ソルバをBareMetal上に再実装して、高速化を測ります。これはいわば、Raspberry Piを1次方程式を解くための専用機として運用するようなものです。このような贅沢な使い方ができるのも、自作スパコンの特権です。

早速、製作を行っていきましょう。

クロスコンパイル開発環境を整える

プログラムのバイナリを動かす環境と、コンパイルを行う環境が異なるコンパイルのことを、クロスコンパイルといいます。 BareMetal開発を行うには専用のクロスコンパイル開発環境が必要なので、これから用意していきます。

以降の作業をRaspberry Piで行うと、最低でもまる1日以上かかってしまいます。なので性能の良い別のマシン上（以降、母艦）にクロスコンパイル環境を用意します。

環境を合わせるため、母艦のOSはDebian Linux(wheezy） x86_64版を用います。他の環境を使われている方は、仮想マシン上にこの環境を用意して行ってください。

Crosstool-NGのインストール

クロスコンパイル環境作成は、Crosstool-NGというツールを使って行います。以下のコマンドを実行してインストールしてください。

$ wget http://crosstool-ng.org/download/crosstool-ng/crosstool-ng-1.20.0.tar.bz2
$ tar jxvf crosstool-ng-1.20.0.tar.bz2
$ cd crosstool-ng-1.20.0
$ ./configure
$ make
$ sudo make install

クロスコンパイラのビルド

ここからはCrosstool-NGを使ってクロスコンパイル環境を作成します。作成に必要な設定ファイルを以下からダウンロードして、母艦のホームディレクトリに保存してください。

hard_dot_config

ダウンロード完了後、以下のコマンドを実行すると、クロスコンパイル環境のビルドが始まります。環境にもよりますが、ビルド完了までは大体３０分以上かかります。気長にお待ちいただければと思います。

$ mkdir -p ~/ct-ng-rpi/src
$ cd ~/ct-ng-rpi
$ ct-ng arm-unknown-eabi
$ mv ~/hard_dot_config .config
$ ct-ng build

パスの追加

作成の終わったクロスコンパイル環境は ~/cross/arm-hardfloat-eabi 以下に保存されます。このツールを使えるよう、システムにパスを通します。以下のコマンドを実行してください。

$ export PATH=$PATH:$HOME/cross/arm-hardfloat-eabi/bin
$ echo 'export PATH=$PATH:$HOME/cross/arm-hardfloat-eabi/bin' >> ~/.bashrc

これで開発環境は整いました。

プログラムのクロスコンパイルとインストール

BareMetal環境上ではOSのもつ資源が一切ないので、初期状態ではprintfどころかC言語自体がまともに動きません。

この問題を解決するには自分で初期化コードを書く必要があります。しかし、この作業は今回の本質でないので省略して、著書のコードをベースに開発を行います (この内容について興味のある方は、著書を読んでいただけると嬉しいです)。

以下に著書のコードをベースとして、今回作成したLU分解の連立一次方程式ソルバをMyLUsolver関数として実装したコードを示します。次のリンクからダウンロードして下さい。

LU_baremetal_normal.zip

クロスコンパイル

先ほどダウンロードしたコードをクロスコンパイルして、 BareMetal Raspberry Piで動くコードを生成します。

コードにはMakefileを用意してあるので、makeコマンドでコンパイルは完了します。

$ make

実際に動かすバイナリの名前は rpi-micon.imgです。

インストール

Raspberry Piは起動時に、第1パーティションの設定ファイル config.txt を読み込んで、各種初期設定を行います。この設定ファイルを書き換えると、起動時に読み込むプログラムを書き換えられます。

ここからはその設定ファイルを書き換え、自作コードが動くように設定を行います。 Raspberry Piに刺していたSDカードを抜き、パソコンに接続してください。

するとこのようにBOOT（第1パーティション）とRECOBERYの2つのパーティションが見えるはずです。

BOOTパーティションを開くと、中にconfig.txtが見つかります。このファイルの末尾に、以下の1行を追記してください。

kernel=rpi-micon.img

これでRaspberry Piは、第1パーティションのrpi-micon.imgを起動するようになります。

最後にrpi-micon.imgを母艦からSDカードに移せばインストールは完了です。

シリアル通信

BareMetal環境上の標準出力はRaspberry PiのUART（シリアル通信）に出ていて、信号レベルは3.3Vです。これを母艦で読み取るため、USB-シリアル変換ケーブルとRaspberry Piを接続します。

※シリアルケーブルはFTDI社のケーブル（秋月電子）がお勧めです。

Raspberry PiのUARTはP1コネクタのGPIO14番からTX、GPIO1510番からRXが出ています。

上記の図と、USB-シリアル変換ケーブルの仕様を参照して
・Raspberry PiのTXをUSB-シリアル変換ケーブルのRX
・Raspberry PiのRXをUSB-シリアル変換ケーブルのTX
・Raspberry PiのGNDをUSB-シリアル変換ケーブルのGND
に接続してください。

シリアル通信を行うプログラムは、Windowsの方はTeraTerm、Macの方はCoolTermをおすすめします。ボーレートは115200bpsに設定してください。

BareMetal版の性能測定

SDカードをRaspberry Piに刺し電源を起動するとプログラムが動作し、連立一次方程式の解と実行時間を出力します。

最初に比較のため、Linux上で動かした連立一次方程式ソルバの結果を示します。コードは以下からダウンロードして下さい。

LU_normal_192x192.zip

$ gcc -O2 -o LU_normal_192x192 LU_normal_192x192.c
$ ./LU_normal_192x192

exec time: 0.380212(sec)

結果は380ミリ秒でした。これに比べてBareMetalでは、はたしてどれだけ性能が向上するでしょうか。以下に結果を示します。

exec time: 458917(us)

459ミリ秒とLinux版に対して20％ほど余計に時間がかかっています。

……忘れていました。 OSが無いということは、ただ開発が不便なだけではないのでした。起動直後のCPUは、ハードウェア浮動小数点演算器（VFP）やキャッシュなど、様々な機能がOFFになっています。そのため起動時にはこれらの機能を適切に設定することで、ハードウェアの性能を存分に発揮できる状態にしてあげなければいけません。それを行うコードもOSのカーネルに組み込まれています。

私は今回、VFPの設定は行いましたが、キャッシュの設定を行っていませんでした。そのためどうやら、キャッシュメモリの設定が行われているLinuxの方は、命令かデータがキャッシュに残っており、その差が実行時間の差に現れたのだと考えられます。

リベンジ——VFPを効率的に用いる最適化

このままでは悔しいので、BareMetal版コードの方をいじって高速化してみます。

Raspberry PiのCPUにはVFP(v2)というハードウェア浮動小数点演算器がついており、これを使うと浮動小数点演算を高速に処理できます。早速やってみましょう。

通常、ARMのCPUで計算を行うには、CPUの持つ汎用レジスタ（r0,r1,...,r12）を用います。しかし、VFPはCPUとは独立した装置なので、VFPの計算はVFPが持つレジスタ間のみで行われます。そのため計算を行うには、値をメモリからVFPのレジスタにコピーし、結果をメモリにストアする作業が行われます。

一般的に、メモリのロード・ストアや、異なる装置間の値交換は、レジスタのみを用いた演算に比べ低速です。しかし、この遅い処理を一度にまとめて行うと、パイプラインが効いて速く処理できる場合があります。この恩恵を最大限に受けられるようにしましょう。

ループのアンローリング

Raspberry PiのVFPは16本の64bitレジスタ（d0~d15）を持っているので、最大16個の値を一度にまとめて処理して効率化できます。これがどこに適応できるかというと、LU分解の行列Lを処理した以下の部分と、

// 行列Lの処理
double div_tmp = 1.0 / A[k][k];

/////////////////////////////////
////// 適応できる場所ここから //////
/////////////////////////////////
for (i=k+1; i<n; i++) {
    // iで回す行方向ループ

    // l_ik = a_ik / a_kk
    A[i][k] = A[i][k] * div_tmp;
}
/////////////////////////////////
////// 適応できる場所ここまで //////
/////////////////////////////////

並列化の時にも触った、行列Uを処理する以下の部分です。

// 行列Uの処理
for (j=k+1; j<n; j++) {
    // jで回す行方向のループ

    // Tips: ここで値を変数に入れたほうが
    // ループ内で参照するよりメモリアクセス回数が減ってお得
    double a_jk = A[j][k];

    /////////////////////////////////
    ////// 適応できる場所ここから //////
    /////////////////////////////////
    for (i=k+1; i<n; i++) {
        // iで回す列方向のループ
        A[j][i] = A[j][i] - A[k][i] * a_jk;
    }
    /////////////////////////////////
    ////// 適応できる場所ここまで //////
    /////////////////////////////////
}

ここを一度にまとめて行うよう、次のように書き換えます。

/ 行列Lの処理
double div_tmp = 1.0 / A[k][k];
for (i=k+1; i<n-12; i+=12) {
    // iで回す行方向ループ

    double d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,d15;

    // l_ik = a_ik / a_kk
    d4 = A[i][k];
    d5 = A[i+1][k];
    d6 = A[i+2][k];
    d7 = A[i+3][k];
    d8 = A[i+4][k];
    d9 = A[i+5][k];
    d10 = A[i+6][k];
    d11 = A[i+7][k];
    d12 = A[i+8][k];
    d13 = A[i+9][k];
    d14 = A[i+10][k];
    d15 = A[i+11][k];

    d4 *= div_tmp;
    d5 *= div_tmp;
    d6 *= div_tmp;
    d7 *= div_tmp;
    d8 *= div_tmp;
    d9 *= div_tmp;
    d10 *= div_tmp;
    d11 *= div_tmp;
    d12 *= div_tmp;
    d13 *= div_tmp;
    d14 *= div_tmp;
    d15 *= div_tmp;

    A[i][k] = d4;
    A[i+1][k] = d5;
    A[i+2][k] = d6;
    A[i+3][k] = d7;
    A[i+4][k] = d8;
    A[i+5][k] = d9;
    A[i+6][k] = d10;
    A[i+7][k] = d11;
    A[i+8][k] = d12;
    A[i+9][k] = d13;
    A[i+10][k] = d14;
    A[i+11][k] = d15;
}
for(;i<n;i++){
    A[i][k] = A[i][k] * div_tmp;
}
// 行列Uの処理
for (j=k+1; j<n; j++) {
    // jで回す行方向のループ

    // Tips: ここで値を変数に入れたほうが
    // ループ内で参照するよりメモリアクセス回数が減ってお得
    double a_jk = A[j][k];

    for (i=k+1; i<n-6; i+=6) {
        // iで回す列方向ループ
        double d4,d5,d6,d7,d8,d9;

        // l_ik = a_ik / a_kk
        d4 = A[k][i+0];
        d5 = A[k][i+1];
        d6 = A[k][i+2];
        d7 = A[k][i+3];
        d8 = A[k][i+4];
        d9 = A[k][i+5];

        d4 *= a_jk;
        d5 *= a_jk;
        d6 *= a_jk;
        d7 *= a_jk;
        d8 *= a_jk;
        d9 *= a_jk;

        A[j][i+0] -= d4;
        A[j][i+1] -= d5;
        A[j][i+2] -= d6;
        A[j][i+3] -= d7;
        A[j][i+4] -= d8;
        A[j][i+5] -= d9;
    }
    for(;i<n;i++){
        A[j][i] = A[j][i] - A[k][i] * a_jk;
    }
}

このようにループを展開して高速化を行う手法を、ループのアンローリングといいます。

行列Lで行っていた計算は、

A[i][k] = A[i][k] * div_tmp;

このようなdiv_tmpと行列Aの要素1つの単純な掛け算なので、使用するVFPのレジスタ16以下に収まるよう、12回分(1+12 < 16)を一度に行うよう展開しています。行列Uの計算と更新は、

A[j][i] = A[j][i] - A[k][i] * a_jk;

このようにa_jkと行列Aの要素を掛けた後、行列Aの別の要素へ引き算を行っています。これはa_jkに１つ、行列Aの要素に2つレジスタを使うので、6回分（1+2*6 < 16）を一気に行うよう展開しています。

ループのアンローリング効果の確認

ループのアンローリングによる効果は、逆アセンブル結果を見れば明らかです。以下にアンローリング後のコードを示すので、ダウンロードしてmakeを行ってください。

LU_baremetal_unrolling.zip

逆アセンブル結果は rpi-micon.disas に出力されます。このファイルを、アンローリング前と後で見比べてみてみます。

行列Lを処理していた部分の、アンローリング前のコードを示します。

809c:   ed937b00    vldr    d7, [r3]
80a0:   ee277b06    vmul.f64    d7, d7, d6
80a4:   ed837b00    vstr    d7, [r3]
80a8:   e2833c06    add r3, r3, #1536   ; 0x600
80ac:   e15c0003    cmp ip, r3
80b0:   1afffff9    bne 809c <myLUsolver+0x34>

このコードはループを使って処理しているので、1変数ずつvldrでロード、vmul.f64で乗算、vstrでストアを繰り返しています。しかしアンローリングを行うと、次のようにロードとストアが一列に並ぶコードになります。

8124:   ed92ab00    vldr    d10, [r2]
8128:   ed919b00    vldr    d9, [r1]
812c:   ed9b8b00    vldr    d8, [fp]
8130:   ed9a0b00    vldr    d0, [sl]
8134:   ed991b00    vldr    d1, [r9]
8138:   ed982b00    vldr    d2, [r8]
813c:   ed953b00    vldr    d3, [r5]
8140:   ed944b00    vldr    d4, [r4]
8144:   ed9c5b00    vldr    d5, [ip]
8148:   ed93bb00    vldr    d11, [r3]
814c:   ed97cb00    vldr    d12, [r7]
8150:   ed906b00    vldr    d6, [r0]
8154:   ee2aab07    vmul.f64    d10, d10, d7
8158:   ee299b07    vmul.f64    d9, d9, d7
815c:   ee288b07    vmul.f64    d8, d8, d7
8160:   ee200b07    vmul.f64    d0, d0, d7
8164:   ee211b07    vmul.f64    d1, d1, d7
8168:   ee222b07    vmul.f64    d2, d2, d7
816c:   ee233b07    vmul.f64    d3, d3, d7
8170:   ee244b07    vmul.f64    d4, d4, d7
8174:   ee255b07    vmul.f64    d5, d5, d7
8178:   ee2bbb07    vmul.f64    d11, d11, d7
817c:   ee2ccb07    vmul.f64    d12, d12, d7
8180:   ee266b07    vmul.f64    d6, d6, d7
8184:   e59d3000    ldr r3, [sp]
8188:   e283300c    add r3, r3, #12
818c:   e58d3000    str r3, [sp]
8190:   e35300b3    cmp r3, #179    ; 0xb3
8194:   e59d3008    ldr r3, [sp, #8]
8198:   ed82ab00    vstr    d10, [r2]
819c:   ed819b00    vstr    d9, [r1]
81a0:   ed8b8b00    vstr    d8, [fp]
81a4:   ed8a0b00    vstr    d0, [sl]
81a8:   ed891b00    vstr    d1, [r9]
81ac:   ed882b00    vstr    d2, [r8]
81b0:   ed853b00    vstr    d3, [r5]
81b4:   ed844b00    vstr    d4, [r4]
81b8:   ed8c5b00    vstr    d5, [ip]
81bc:   ed83bb00    vstr    d11, [r3]
81c0:   e2833b12    add r3, r3, #18432  ; 0x4800
81c4:   ed87cb00    vstr    d12, [r7]
81c8:   e2822b12    add r2, r2, #18432  ; 0x4800
81cc:   ed806b00    vstr    d6, [r0]
81d0:   e2811b12    add r1, r1, #18432  ; 0x4800
81d4:   e28bbb12    add fp, fp, #18432  ; 0x4800
81d8:   e28aab12    add sl, sl, #18432  ; 0x4800
81dc:   e2899b12    add r9, r9, #18432  ; 0x4800
81e0:   e2888b12    add r8, r8, #18432  ; 0x4800
81e4:   e2855b12    add r5, r5, #18432  ; 0x4800
81e8:   e2844b12    add r4, r4, #18432  ; 0x4800
81ec:   e28ccb12    add ip, ip, #18432  ; 0x4800
81f0:   e58d3008    str r3, [sp, #8]
81f4:   e2877b12    add r7, r7, #18432  ; 0x4800
81f8:   e2800b12    add r0, r0, #18432  ; 0x4800
81fc:   daffffc8    ble 8124 <myLUsolver+0xbc>

美しいですね。

連続してロードストアや演算を行っているので、パイプラインが効きやすいコードになっているはずです。

最適化結果の確認

ではこのプログラムを実行して、どれだけ最適化できたか確認してみましょう。以下に結果を示します。

exec time: 417205(us)

アンローリングすると、その分コードサイズは増えてしまいますが、命令の実効効率が良くなります。その証拠に、アンローリング後のプログラムの実行時間は約40msも短くなりました。

しかし、Linux上で動かした結果には後一歩及びませんでした。リベンジ失敗です。

BareMetal版が勝つためには……

BareMetal版がLinuxで動かした結果に勝つためには、すくなくともLinuxで行っているハードウェアの設定を、BareMetalでも行う必要があるでしょう。しかし、その設定を行うため、CPUの資料やコードを読むコストを考えると、BareMetalでやるより、Linuxカーネルの設定を修正したほうが速いでしょう。

BareMetal版が優っている点としては、電源投入後すぐに実行され、結果が出ることでしょうか。なので、電源投入から解が出るまでの時間で競えば、BareMetal版の勝利です。

今回のまとめ

・並列プログラミングの勉強には連立一次方程式ソルバがおすすめ
・並列化しても高速にならないどころか、逆に遅くなることもある
・OSの資源は偉大
・ループのアンローリングを適切に行うと、処理を高速化できる

次回は通信コストを削減して、再度、並列処理に挑戦する予定です

著者プロフィール

西永俊文　1992年生まれの情報系大学生。一人前のエンジニアを目指して日々勉強中。PDA全盛期の影響から、組み込み機器が好き。今欲しい物は、現代の技術で作られた物理キーボード付きで4インチ以下の小さなスマートフォン。著書に『BareMetalで遊ぶ Raspberry Pi』がある。
〈四コママンガ：囚〉

OLYMPUS DIGITAL CAMERA — 著者近影（手にしているのは今回製作したRasPiクラスター初号機）

連載目次

2014-12-01 RaspberryPiクラスタ製作記第0回「スパコンって本当にスゴイの？」
2014-12-01 PHASE/0のインストール方法
2014-12-16 RaspberryPiクラスタ製作記第1回「スパコンを作ろう」
2014-12-16 Raspberry Piの初期設定
2014-12-29 RaspberryPiクラスタ製作記第2回「スパコンで遊ぼう」
2014-12-29 HPLのインストール方法
2015-02-11 RaspberryPiクラスタ製作記第3回「並列プログラミング」
2015-02-11 MPIを用いた並列プログラミングの概要
2015-02-11 LU分解アルゴリズムのおさらい
2015-03-01 RaspberryPiクラスタ製作記第4回「BareMetal並列計算」

HPLのインストール方法

Raspberry PiへのHPLインストールの手順を説明します。

1. 依存ライブラリのインストール

HPLのコンパイルに必要なライブラリをインストールします。 Raspberry Pi上で以下のコマンドを実行してください。

$ sudo aptitude -y install libblas3 libblas-dev libatlas libatlas-base-dev libc6-dev

2. HPLのダウンロードと解凍

HPL配下のページで配布されています。
http://www.netlib.org/benchmark/hpl/

今回は2014年12月15日現在最新の hpl-2.1 を用いて計算を行います。以下のコマンドでファイルをダウンロード、解凍してください。

$ wget http://www.netlib.org/benchmark/hpl/hpl-2.1.tar.gz
$ tar zxvf hpl-2.1.tar.gz

解凍して出てきたディレクトリは、次に行うmake作業を楽にするため、ホームディレクトリに移動します。以下のコマンドを実行してください。

$ mv hpl-2.1 ~/hpl

3. Makefileの設定

HPLをmakeするため、ライブラリやコンパイラの場所をMakefileに設定します。本来の手順では以下のコマンドのように ./setup ディレクトリのテンプレートをhplのルートディレクトリにコピーし、設定を行っていきます。

$ cp setup/Make.Linux_ATHLON_CBLAS ./Make.Linux_ARM_CBLAS

しかしこの設定の手順は非常に大変なので、今回は私が設定したMakefileを配布し、これを用いてmakeを行います。以下のURLよりファイルを取得してください。

4. make

HPLのmakeは、先ほどダウンロードしたMakefileの拡張子を引数archに渡して行います。以下のコマンドを実行してください。

$ make arch=Linux_ARM_CBLAS

5. コピー

これまでの手順をすべてのRaspberry Pi上で行うのは大変です。そこで、すでにコンパイルの終わったものを、rsyncで他のRaspberry Piにコピーします。

以下にrpicluster02をコピー先としたコマンドを示すので、このコマンドの”rpicluster02″の部分を各マシンごとに変更しつつ、すべてのマシン宛に行ってください。

$ rsync -ah -e ssh ~/hpl pi@rpicluster02:~/

これでインストールは完了です。

RaspberryPiクラスタ製作記第2回「スパコンで遊ぼう」

第２回　Raspberry Piスパコンで遊ぼう

自作したばかりのスーパコンピュータ上で、第0回で動かした第一原理電子状態計算ソフトウェアPHASE/0を動かして遊んでみましょう。

この遊びの目的は以下の2つです。

・自作スパコンでも「京」で動いている計算プログラムが動くことを確認する
・並列化により処理がどれだけ高速化されるかを確認する

「並列化により処理がどれだけ高速化されるか」に対して、まずあえてもっとも単純に予想すると「並列計算ノード数を増やすと、それに比例して性能も上昇する」となります。ノード数を倍にするということは計算に使えるリソースが倍になるということなので、ノード数に比例して処理速度も上がっていくのではないでしょうか。この予想が正しいかどうかを、これから検証していきます。

PHASE/0のインストール

今回使用するRaspberry PiはARMアーキテクチャのCPUを搭載しているため、第0回でx86マシン用にコンパイルしたバイナリを実行することができません。そこで再度ソースコードからのコンパイルを行います。

最初にPAHSE/0のコンパイルと動作に必要なパッケージをインストールします。以下のコマンドを実行してください。

$ sudo aptitude -y install gfortran mpich2 libc6-dev

PHASE/0のコンパイルは基本的に ./install.sh を実行するだけですが、こちらで試したところ、今回使用したバージョンのインストールスクリプトはRaspbianでは正常に動作しませんでした。そこで、このインストールスクリプトが正常動作する環境（MacOSX)で動作させてコンパイルに必要なMakefileを生成し、これをRaspberry Pi上に持ち込むことでコンパイルを可能としました。このMakefileを以下のURLより取得してください。

Makefile

取得したMakefileを phase0_2014.03/src_phase/ 以下に置いて、そのディレクトリ上で make コマンドを実行するとコンパイルが開始されます。エラーなくコンパイルが終了すれば、PHASE/0のインストールは完了です。

計算打ち切り時間の書き換え

PHASE/0の入力ファイル中には、計算をどれだけで打ち切るか（cpumax）の設定が書かれています。今回使用するRaspberry Piは第0回で使用したマシンに比べて非常に性能が低いため、デフォルトの打ち切り時間では計算が終了しません。そのためこの設定を修正して、打ち切り時間を延長します。

samples/Si8/input_scf_Si8.data を開いて、以下の3600secとなっている部分を1dayに修正してください。

Control{
		# 修正前
        cpumax = 3600 sec     ! {sec|min|hour|day}
		# 修正後
        cpumax = 1 day     ! {sec|min|hour|day}
        condition=initial
}

MPIノードの設定

PHASE/0のプログラムはMPI方式でデータをやり取りする、mpichというアプリケーションを用いて並列計算を行います。このmpichは、自分のマシンと一緒に計算を行うマシン（ノード）のホストネーム（またはIPアドレス）を設定ファイルに記述し、それを実行時に読み込ませることで並列計算が可能となります。以下にその設定ファイルを記すので ~/mpdhosts に保存してください。

rpicluster01
rpicluster02
rpicluster03
rpicluster04
rpicluster05
rpicluster06
rpicluster07
rpicluster08
rpicluster09
rpicluster10
rpicluster11
rpicluster12
rpicluster13
rpicluster14
rpicluster15
rpicluster16

処理時間の計測

インストールの終わったPHASE/0を使い、ノード数を1,4,9,16と変化させて、その処理時間を計測します。

1ノード1並列

最初は動作確認のため１並列で動かしてみましょう。以下のコマンドで計算を開始してください。

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=1 ../../bin/phase -ne=1 -nk=1
real    32m59.518s
user    29m53.770s
sys     3m2.960s

実行時間は実時間で約33分と、第0回と比べて非常に時間のかかる結果となりました。今回使用したRaspberry Piはクロックが700MHzの32bitマシンなので、妥当なスコアでしょう。

4ノード4並列

次にRaspberry Piクラスタマシンのノードを4台使用し、計算を4並列で回してみます。予想のとおりであれば、CPUの数が4倍になるので実行時間も1/4になるはずですが、果たしてどうなるでしょうか。

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=4 ../../bin/phase -ne=2 -nk=2
real    16m14.437s
user    13m50.940s
sys     2m20.780s

実行時間は実時間で約16分、1並列の約半分となりました。 1/4とまでは行きませんでしたが、1/2でもまずまずのスコアでしょう。

この調子で9,16並列で計算を行い、その性能の上昇を見ていきます。

9ノード9並列

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=9 ../../bin/phase -ne=3 -nk=3
real    11m22.431s
user    9m41.020s
sys     1m39.270s

16ノード16並列

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=16 ../../bin/phase -ne=4 -nk=4
real    9m25.834s
user    7m11.440s
sys     2m11.720s

計算結果の考察

1,4,9,16並列の計算時間の結果をまとめたグラフを示します。

予想ではノード数を上げると計算に使えるCPUの数が増えるので、ノード数に比例して実行時間は減っていくはずでした。しかし実際は、4並列の時の実行時間は一気に約1/2となったものの、その後の9並列では1並列に比べ約1/3と伸びが悪くなり、16並列ではより一層伸びが悪くなっています。

一体どうして、計算時間はノード数に比例して下がっていかなかったのでしょうか。

アムダールの法則

計算時間がノード数に比例して下がらない理由ついて、第1回と同じように夏休みの宿題にたとえて説明します。第1回では、夏休みの宿題として計算問題を100題だされたとしても、 100人の協力者が1問づつ解いて最後に答えを共有すれば、一瞬で解くことができると言っていました。この計算問題を100人で解く手順を改めて示すと、以下の5手順に分けられます。

1.計算を100個に分割する
2.分けた計算を担当する協力者に解いてもらうよう指示を出す
3.それぞれ協力者が計算する
4.結果を集める
5.結果を問題集に書き込む

こうしてみると、100人の協力者が同時に作業を行える部分、つまり並列処理できる部分は手順3.だけです。それ以外の部分は指示を出す人が1人で行う必要がある部分、つまり並列処理することができない部分です。

この例からわかるように、処理の中には並列に処理することにより高速化できる部分と、そうでない部分が存在します。そのためノード数を上げることによる処理時間の減少は、一定のラインを超えると並列処理により高速化できない部分がネックとなり、処理時間の減少の伸びが非常に悪くなります。これが今回私たちの実験で起こった現象です。

このノード数（プロセッサ）と処理速度の関係には「アムダールの法則」という法則があります。アムダールの法則は、並列化できない部分の割合をs、プロセッサ数をnとすると、次の式で処理速度の向上比を計算することができます。

実際にこの式を使い、s=0.001、n=1~30000で計算した結果をグラフに示します。

性能評価とコストの比較

コンピュータの性能には、理論性能と実行性能の2つがあります。

理論性能　ハードウェアの性能から単純に計算した最大性能。例えば700MHzで動作するRaspberry Piは、1クロックで1回浮動小数点演算が行えるので、理論性能700MFLOPSとなる。

実効性能　何らかのベンチマークソフトを動かして得た性能。

理論性能はあくまで理論的な最大性能なので、実効性能とは離れた結果となります。先ほど私たちの自作スパコン上でPHASE0を動かした結果からもわかると思います。しかし、具体的にどれだけの性能を持っているかは、まだ調べていませんでした。

性能の他にもうひとつ気になるのは、そのお値段です。現代はコストパフォーマンスを重視する傾向にあります。そのため、どれだけスパコンの性能が良くとも製作費やランニングコストが高すぎては実用的ではありません。

ところで、Raspberry PiのCPUアーキテクチャはARMというものでした。このARMのCPUは性能に対して消費電力が低いという特徴があり、身近なところでは携帯ゲーム機やスマートフォンで使用されています。そのためもしかすると、消費電力――つまりランニングコストの面では、スパコン「京」に勝つことができるのではないでしょうか。

ここからはRaspberry Piスパコンの実効性能とそのコストを調べ、最後にスパコン「京」と比較してみましょう。

実効性能計測

Linpackは、コンピュータの実効性能の計測（ベンチマーク）を行うプログラムのひとつです。 TOP500ではこれを並列化した「HPL」(High-Performance Linpack)を用いて、ベンチマークスコアを測っています。私たちのRaspberry Piクラスタマシンも、このHPLを用いて実効性能を計ってみましょう。

HPLのインストールについては、別のページで説明します。

HPLは処理する問題サイズ、ブロックサイズ等のパラメータ変更することで最適化を行い、スコアを上げることができます。しかし最適なパラメータを探すのは非常に大変なので、今回は私の試した中で一番よいスコアの得られた設定を利用します。

以下のファイルをダウンロードし、~/hpl/Linux_ARM_CBLAS/以下へ上書きコピーしてください。

HPL.dat

HPLの実行

実行方法はPHASE0とほとんどおなじですが、こちらはノードの指定にホストネームを使うとエラーが発生するので、 IPアドレスを使ってノードを指定します。本環境でのrpicluster01-16のIPアドレスは 192.168.240.101-116なので、16ノード計算を行うため以下のコマンドを実行します。

$ cd ~/hpl/bin/Linux_ARM_CBLAS
$ mpirun -np 16 -host 192.168.240.101,192.168.240.102, \
192.168.240.103,192.168.240.104, \
192.168.240.105,192.168.240.106, \
192.168.240.107,192.168.240.108, \
192.168.240.109,192.168.240.110, \
192.168.240.111,192.168.240.112, \
192.168.240.113,192.168.240.114, \
192.168.240.115,192.168.240.116  \
./xhpl

ベンチマーク結果

=================================================================
T/V                N    NB     P     Q               Time                 Gflops
-----------------------------------------------------------------
WR11C2R4       25600    96     4     4            4211.62              2.656e+00

HPLの結果は2.656GFLOPSとなりました。このスコアは、今から約30年前の1987年に発表されたスパコン「HITAC S-810」（最大3GFLOPS）に匹敵するスコアです。
　参考URL: http://museum.ipsj.or.jp/computer/super/index.html
当時は巨大なシステムで実現していた性能が、今ではこんな小さなボード16枚で得られてしまいます。私は当時まだ生まれていませんが、技術の進歩が肌で感じられました。

ここからは自作Raspberry Piクラスタマシンの制作費とランニングコストを調べていきます。

製作費は11万円

今回製作に使用した部品とそのお値段の一覧を示します。なお、製作の人件費はないものとします。

部品	数	単価(円)
Raspbebrry Pi TypeB+	16	3,900
micro SDカード 8GB	16	800
USB-MicroBケーブル	16	600
2ポートUSB充電器	8	800
LANケーブル	16	350
24ポートスイッチングハブ	1	7,800
9ポート電源タップ	2	1,000
M2.6連結スペーサーネジ	72	10
小型メタルラック	1	600
結束バンド	1	100
ファンシー金網	1	100
金網用フック付き小物いれ	1	100

以上を合計すると、製作にかかったお金の合計は約11万円でした。

電気代（ランニングコスト）

USB充電器（Raspberry Piの外部電源）とUSBケーブルの間に簡易テスターを仕込んで、Raspberry Piクラスタマシンの消費電力を調べます。

結果は１台あたり200mA。Raspberry Piには負荷に応じてクロックを変えるような機能は無いため、ピーク時もこの値です。よって１台あたりの消費電力は 5(V) x 0.2(A) = 1(W)、16台でも16Wです。その他、スイッチングハブやUSB充電器のロスを含めると、このクラスタマシンの総消費電力は大体20Wとなります。

この消費電力から電気料金を計算すると、1日あたり約10円となります。製作費同様、メンテナンスにかかる人件費は0円とすると、ランニングコストは1日あたり10円、1年動かし続けても3560円です。

京との比較

比較のため、スパコン「京」のコストについても調べてみましょう。京のコストについての資料としては、平成25年1月に文部科学省からの資料があります。この資料によると、京の制作費は合計793億円（施設建設費用は含めず）、運用費は平成24年のデータで1年あたり97億円です。

ここまで調べたRaspberry Piクラスタと京の制作費を表に示します。

	実効性能	制作費（円）	運用費（円/年）
京	10PFLOPS	793億	97億
RasPiクラスタ	2.656GFLOPS	11万	3560

これだけではわかりづらいので、性能をコスト（製作費＋運用費）で割った、1FLOPSあたりのコストを計算して比較します。

もしかするとコストの面では勝てるのではと考えていましたが、仮に100年運用したとしても京を追い抜くことはできませんでした。技術の粋を集めて作られた京は、個人のレベルでは到底届かないものであることを再認識することが出来ました。

まとめ

今回はRaspberry Piスパコン制作のソフトウェアの整備と、その上でPAHSE0の動作を行い、ノード数と計算時間の関係について調べてみました。この内容を以下に簡単にまとめます。

・ノード数増加による処理速度向上は一定のところで伸び悩む
・アムダールの法則はプロセッサ数と性能向上比の関係を教えてくれる
・京のコストパフォーマンスは自作Raspberry Piスパコンと比べると圧倒的に良い

著者プロフィール

西永俊文　1992年生まれの情報系大学生。一人前のエンジニアを目指して日々勉強中。PDA全盛期の影響から、組み込み機器が好き。今欲しい物は、現代の技術で作られた物理キーボード付きで4インチ以下の小さなスマートフォン。著書に『BareMetalで遊ぶ Raspberry Pi』がある。

連載目次

↑前のページ／次のページ↓

LU分解アルゴリズムのおさらい

LU分解を使って連立一次方程式を解く方法

LU分解の方法

MPIを用いた並列プログラミングの概要

並列計算機の種類とMPI

MPI関数の概要

MPIを用いた並列HelloWorld

MPI関数を用いたデータ通信サンプルプログラム

RaspberryPiクラスタ製作記 第3回「並列プログラミング」

第３回 並列計算プログラムをつくろう

連立一次方程式をプログラムで解く

行列の定義

LU分解のコードを作る

完成

LU分解コードを並列化しよう

性能測定

BareMetalで作る連立一次方程式ソルバ専用機

クロスコンパイル開発環境を整える

プログラムのクロスコンパイルとインストール

BareMetal版の性能測定

リベンジ——VFPを効率的に用いる最適化

ループのアンローリング効果の確認

最適化結果の確認

BareMetal版が勝つためには……

今回のまとめ

著者プロフィール

連載目次

HPLのインストール方法

1. 依存ライブラリのインストール

2. HPLのダウンロードと解凍

3. Makefileの設定

4. make

5. コピー

RaspberryPiクラスタ製作記 第2回「スパコンで遊ぼう」

第２回 Raspberry Piスパコンで遊ぼう

PHASE/0のインストール

計算打ち切り時間の書き換え

MPIノードの設定

処理時間の計測

計算結果の考察

アムダールの法則

性能評価とコストの比較

実効性能計測

ベンチマーク結果

製作費は11万円

電気代（ランニングコスト）

京との比較

まとめ

著者プロフィール

連載目次

RaspberryPiクラスタ製作記第3回「並列プログラミング」

第３回　並列計算プログラムをつくろう

RaspberryPiクラスタ製作記第2回「スパコンで遊ぼう」

第２回　Raspberry Piスパコンで遊ぼう