数値計算レポート

G99P102-8 永田貴彦

出題日 :2001 年 6 月 6 日提出期限:2001 年 6 月 20 日提出日 :2001 年 6 月 20 日

目次

1 問題 2

2 環境の構築 3 2.1 CLAPACK for Win32(version2.0) の導入 . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 CLAPACK version 3.0 の導入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 標準 BLAS の導入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 最適化 BLAS の導入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 問題の考察 13 3.1 問題 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2 問題 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

4 感想 23

1 問題

問題 1 行列の積のプログラムを使って、 dgemm (...) を作れ。 (教科書「Linux 数値計算ツール」 P.137 参照) CLAPACK をインストールしながら仕様をみよ。

問題 2 Ax = b を解くための LAPACK の関数 dgesv (...) を作れ。 (教科書「Linux 数値計算ツール」 P.147 参照) また、 BLAS を高速化したら高速になる形で実装せよ。つまり BLAS の関数を使用して作る。

2 環境の構築

2.1 CLAPACK for Win32(version2.0) の導入

CPU  AMD Athlon 900MHz(FSB100)
Memory  PC133 SDRAM 384MB
Chipset Asus A7V133-C (VIA KT133A + VT82C686B)
OS  Windows2000 Professional + Service Pack1
Compiler  Microsoft Visual C++ 6.0

表 1: プログラム作成・実行環境

まずは、この表の環境に LAPACK をインストールすることから始めます。本当ならば、教科書にあわせて Linux を使ったほうが分かりやすいとは思いましたが、残念ながら Linux を持ってませんので、 Windows の環境で行うことにしました。今回利用する LAPACK は CLAPACK という C で書かれたものであるため、 Windows でもうまく動いてくれるだろうと思ったからです。

最初に、 Netlib の CLAPACK のページ^{1} から、 CLAPACK をダウンロードすることにしました。この説明には、 version 3.0 であると書かれており、すぐ下に、 version 2.0 の CLAPACK もあったのですが、バージョンが新しいほうが一般的にはいいので、 version3.0 の方をダウンロードすることにしました。また、 CLAPACK for Win32 というなにやら Windows 用のようなものもありましたが、 version 2.0 のものしかありませんでした。しかし、 Windows を使っていることもあり、これもダウンロードしておくことにしました。したがって、ダウンロードしたものは、 CLAPACK version 3.0 の clapack.tgz と CLAPACK for Win32 の CLAw32.zip です。

次に、ダウンロードしてきた CLAw32.zip を解凍して、中を確認したところ、 Release というディレクトリに、 clapack.lib というすでにコンパイル済みのものが入っていました。これならわざわざコンパイルしなくてもよさそうです。

これでとりあえずコンパイルは完了したので、今回の課題に使うソースファイルを探してみました。探すものは、 dgemm と dgesv の 2 つです。検索したところありませんでした。どうやらこの Windows 用のパッケージにはソースファイルが入っていないようです。したがって、 clapack.tgz の方を探してみることにしました。その結果、 dgemm は、 CLAPACK/BLAS/SRC/dgemm.c にありました。この関数は、教科書にも、 LEVEL3 BLAS のルーチンとして紹介されているから、この C 言語ソースファイルに間違いないでしょう。もう一方の dgesv は、 CLAPACK/SRC/dgesv.c にありました。この関数はどうやら LAPACK の関数のようです。したがって、 BLAS を利用してこの関数を作れという課題の意味が分かりました。きっと、このソースファイルを参考にして作り、 BLAS を利用したプログラムの書き方に慣れてもらおうという意図でしょう。しかし、丸々コピーして終わりそうな気もします。

どうやら、今回のプログラムでは、 LAPACK を全てインストールする必要はないようです。 BLAS さえ C のプログラムから参照できるようになってればよさそうです。

まず、現在ある標準 BLAS で LAPACK を使えるようにしてから、最適化 BLAS でも使えるようにしてみましょう。先ほどのコンパイル済みのものをどのように使うかいまいちよく分からなかった^{2}ので、説明書通りに自分でコンパイルしながら、理解していくことにします。インストールの手順を読むと、このキット (CLAw32.zip) にはソースコードが含まれていないので、 clapack2.tgz(CLAPACK version2.0) を解凍した上に、上書きしろと書かれていました。どうやら、上書きすることで Visual C++ 6.0 でコンパイルできるようになるらしいです。ただ、 f2c.h というファイルを、それぞれのフォルダに入れてあげないと、エラーが沢山でてきました。それぞれのフォルダにも f2c.h というファイルがありますが、内容が書かれていないためらしいです。 VC++用のプロジェクトファイルがあるので、それを使ってコンパイルが簡単にできます。どうやら、オブジェクトファイルを自動的にライブラリファイルにまとめてくれるらしいです。

1http://www.netlib.org/clapack/

2どうやら、 CLAPACK version2.0 のライブラリと version3.0 のヘッダーファイルを組み合わせていたため、おかしかったらしい。これをあきらめるまで 2〜3 時間も費やしてしまった…

説明書どおり、 libF77→libI77→BLAS→CLAPACK の順番にコンパイルしていったところ、 libF77.lib、 libI77.lib、 blas.lib、 clapack.lib ができました。これらのライブラリはまとめて CLA-PACK/lib/に置き、 VC++のメニューから、「ツール → オプション → ディレクトリ → ライブラリファイル」を選択し、そこのディレクトリへのパスを追加しました。

動作を確認するため、教科書 P.147 のプログラムを実行してみました。 VC++において先ほどのライブラリとリンクさせるためには、メニューの「プロジェクト → 設定 → リンク」を選び、そこにある「オブジェクト/ライブラリモジュール」先ほどの 4 つを追加します。本当に 4 つ必要なのかは知りませんが、全部入れとけば問題ないでしょう。次のような計算をさせることになります。

( \begin{matrix}1 & 1 & 1 \\ 3 & 1 & −3 \\ 1 & −2 & −5 \end{matrix} ) x = ( \begin{matrix}1 \\ 5 \\ 10 \end{matrix} )

答えは、次のようになります。

x = ( \begin{matrix}6 \\ −7 \\ 2 \end{matrix} )

前回作成した、 LU 分解のプログラムとも結果を比較してみましたが、同じ結果が出ていました。 CLAPACK のインストールは成功しているようです。

2.2 CLAPACK version 3.0 の導入

しかし、どうせなら CLAPACK version3.0 を使ってみたいものです。そこでこちらは Windows 用はないものの、何とかコンパイルしてみます。 VC++においてライブラリファイルは、メニューの「ファイル → 新規作成 →Win32 Static Library」から作ることができます。試しに、 libF 77.lib から作成してみます。新規作成したあと、「プロジェクト → プロジェクトに追加 → ファイル」を選び、 CLAPACK/F2CLIBS/libF77/に入っているファイルを全て追加しました。その後、コンパイルさせてみたところ、多くのエラーが出てきました。そこで、先ほどコンパイルした Windows 用のものとどこが違うか調べたところ、 f2c.h に次のようなものを加える必要があることが分かりました。

#ifdef _WIN32
#define _COMPLEX_DEFINED
typedef __int64 longint;
// warning C4244: 'xoperx' : conversion from 'xxx ' to 'xxx ',
//                           possible loss of data
// warning C4101: 'xxx' : unreferenced local variable
#pragma warning( disable: 4244 4101 )
#endif

また、 C++では、 f2ch.add のファイルの内容も加える必要があるようです。しかし、まだまだエラーは無くなりません。次は、「インクルードファイルがオープンできません。 ’sys/times.h’: No such file or directory」といわれました。まあ Unix 依存のライブラリなので当たり前でしょう。まずは、定義を次のように変更します。

#ifndef _WIN32
#include "sys/times.h"
#endif

そして、時間を計算する部分を次のようにしました。

#if defined _WIN32
    clock_t t;
    static clock_t t0;
    t = clock();
    tarray[0] = (t - t0) / CLOCKS_PER_SEC;
    tarray[1] = (t - t0) / CLOCKS_PER_SEC;
    t0 = t;
    return tarray[0] + tarray[1];
#else
    struct tms t;
    static struct tms t0;
    times(&t);
    tarray[0] = (t.tms_utime - t0.tms_utime) / Hz;
    tarray[1] = (t.tms_stime - t0.tms_stime) / Hz;
    t0 = t;
    return tarray[0] + tarray[1];
#endif

これは、 Windows 用のものを参考にして、自分で手直ししたので、あっているかは分かりません。ただ、似たようなコードだったので多分あってるでしょう。

次に、 unsigned long long u; というところでエラーになりました。このような宣言はできたのでしょうか？確かにコメントに機種依存って書いてありますが、機種依存はなるべく避けて欲しいものです。仕方ないので、

#ifndef _WIN32
    unsigned long long u;   /* system-dependent */
#else
    unsigned long u;
#endif

にしました。

また次に、「warning C4723: potential divide by 0 return x == 0 ? 1/x : 0;」というエラーがでてきました。確実に 0 で割っています。これも次のようにしました。

#ifndef _WIN32
    return x == 0 ? 1/x : 0;
#else
    return x == 0 ? 0 : 1/x;
#endif

まだまだエラーがでてきます。次は、「error C2065: ’ulongint’ : 定義されていない識別子です。」がでてきました。確かに、 f2c.h のところには、その定義は存在していますが、＃ if 0 で囲まれています。要するにコメントアウトされているわけです。その定義は、 unsigned long long という定義だったので、使う変数を unsigned long に変えました。

#ifndef _WIN32
typedef unsigned long long ulongint;
#else
typedef unsigned long ulongint;
#endif

ここまで変更したらなんとかコンパイルが通り、 libF77.lib ができました。同様にして、 libI77.lib、 blas.lib、 clapack.lib を作成しました。

では、先ほど試した教科書の例を動作させてみたところ、正確に答えがでているようです。しかし、一体どこが変わったのかは定かではありませんが、バージョンが新しくなって気分的にはいいです。これで、 CLAPACK の導入が終わりました。

ここで変更したプログラムの差分や、追加した VC++用のプロジェクトファイルは、 Web 上^{3} にアップロードしてあります。ただ、変更したファイルを配布してよいのか分からないことと、動作確認を少しはしましたがバグがある可能性があるため、 2 週間後には削除します。また、トップページからリンクによってたどることはできません。

2.3 標準 BLAS の導入

では次に、 C のプログラムから BLAS を使ってみます。試してみるものは問題にもなっている、 dgemm (...) です。教科書には、引数が載っていないので、ソースを調べることにします。ソースを確認したところ、次のようになっていました。

int dgemm_(char *transa, char *transb, integer *m, integer *
    n, integer *k, doublereal *alpha, doublereal *a, integer *lda,
    doublereal *b, integer *ldb, doublereal *beta, doublereal *c__,
    integer *ldc)

ぱっと見、すべて参照渡しなので、値を渡すのか、配列を渡すのかは分かりません。また、変数の宣言も特殊なものになっています。これらはどうやら、 f2c.h に宣言されているらしいです。 integer = long int(VC++では int 型) 、 doublereal = double となっています。では、引数の意味を確認します。これらはソースの中に英語で書かれていました。

• dgemm 行列-行列演算を行います。それは次のようなものです。

C = α × op(A) × op(B) + β × C

3http://www.geocities.co.jp/SiliconValley-SanJose/8084/clapack/clapack sub.zip

α, β はスカラーで、 A, B, C は行列です。 op(A) は m×k 行列、 op(B) は k ×n 行列、 C は m×n 行列です。

• transa char[1] 型 TRANSA は、次の規則にしたがって行列積演算に使われる op(A) の形状を指定します。 TRANSA = ’N’ or ’n’, op( A ) = A TRANSA = ’T’ or ’t’, op( A ) = A^{t} TRANSA = ’C’ or ’c’, op( A ) = A^{t} この引数は、関数中で変化しません。 (新しい値に変更はされません)

• transb char[1] 型 TRANSB は、次の規則にしたがって行列積演算に使われる op(B) の形状を指定します。 TRANSB = ’N’ or ’n’, op( B ) = B TRANSB = ’T’ or ’t’, op( B ) = B^{t} TRANSB = ’C’ or ’c’, op( B ) = B^{t} この引数は、関数中で変化しません。

• M int 型 M は、 op(A) 行列と C 行列の row(行) の値を指定します。 M は 0 以上でなければなりません。この引数は、関数中で変化しません。

• N int 型 N は、 op(B) 行列と C 行列の column(列) の値を指定します。 N は 0 以上でなければなりません。この引数は、関数中で変化しません。

• K int 型 K は、 op(A) 行列の column(列) と op(B) 行列の row(行) の値を指定します。 K は 0 以上でなければなりません。この引数は、関数中で変化しません。

• ALPHA double 型 ALPHA は、 α の値を指定します。この引数は、関数中で変化しません。

• A double 型配列 (LDA,ka) 次元ここで、 TRANSA = ’N’ or ’n’ のときは、 ka = k であり、それ以外は、 ka = m です。 TRANSA = ’N’ or ’n’ の場合には、 A 配列の一部分の (m,k) は、行列 A に含まれている必要があります。それ以外の場合は、 A 配列の一部分の (k,m) は、行列 A に含まれている必要があります。この引数は、関数中で変化しません。

• LDA int 型 LDA は、 (サブ) プログラムを呼ぶときに明示する、 A の最初の次元を指定します。 TRANSA = ’N’ or ’n’ のとき、 LDA は少なくとも、 max(1,m) でなくてはいけません。それ以外は LDA は少なくとも、 max(1,k) でなくてはいけません。この引数は、関数中で変化しません。

• B double 型配列 (LDB,kb) 次元ここで、 TRANSA = ’N’ or ’n’ のときは、 kb = n であり、それ以外は、 ka = k です。 TRANSB = ’N’ or ’n’ の場合には、 B 配列の一部分の (k,n) は、行列 B に含まれている必要があります。それ以外の場合は、 B 配列の一部分の (n,k) は、行列 B に含まれている必要があります。この引数は、関数中で変化しません。 • LDB int 型

LDB は、 (サブ) プログラムを呼ぶときに明示する、 B の最初の次元を指定します。 TRANSB = ’N’ or ’n’ のとき、 LDB は少なくとも、 max(1,k) でなくてはいけません。それ以外は LDB は少なくとも、 max(1,n) でなくてはいけません。この引数は、関数中で変化しません。

• BETA double 型 BETA は、 β の値を指定します。もし BETA が 0 であるなら、 C は値を代入しておく必要はありません。この引数は、関数中で変化しません。

• C double 型配列 (LDC,n) 次元 C 配列の一部分の (m,n) は、行列 C に含まれている必要があります。ただし、 BETA が 0 である場合を除きます。その場合は C に値を入れておく必要はありません。この引数は、終了時に α × op(A) × op(B) + β × C の結果の m × n 行列として上書きされます。

• LDC int 型 LDC は、 (サブ) プログラムを呼ぶときに明示する、 C の最初の次元を指定します。 LDC は少なくとも、 max(1,m) でなくてはいけません。この引数は、関数中で変化しません。

• この関数は、 LEVEL3 BLAS のルーチンです。

では、この BLAS 関数を組み込んだ行列積を計算する関数を作ってみます。行列は次のように定義されていました。

struct matrix {
double **data;
int row;
int column;
};

この行列構造体を、 BLAS での引数に変換する関数をまず作りました。

/* matrix から Blas 用の配列に変換 (m,n) 行列 */
bool matrix2blas(struct matrix* mt, double* a,int m, int n) {
int i,j;
if (m <= 0 || n <= 0) return false;
if ((mt->row < m) || (mt->column < n)) return false;
for (i=0;i<m;i++) {
   for (j=0;j<n;j++) {
       a[i+j*m] = mt->data[i][j];
   }
}
return true;
}
/*Blas 用の配列から matrix に変換 (m,n) 行列 */
bool blas2matrix(double* a,int m, int n, struct matrix* mt) {
int i,j;
if (m <= 0 || n <= 0) return false;

if ((mt->row < m) || (mt->column < n)) return false;
for (i=0;i<m;i++) {
    for (j=0;j<n;j++) {
        mt->data[i][j] = a[i+j*m];
    }
}
return true;
}

そして、行列積の演算は次のように書けます。

bool multi_matrix_dgemm(struct matrix* ma,struct matrix* mb,struct matrix* mc) {
char transa[1];
char transb[1];
int m,n,k;
double alpha,beta;
double* a;
double* b;
double* c;
int lda,ldb,ldc;
if (!is_matrix(*ma) || !is_matrix(*mb)) return false;
if (ma->column != mb->row) return false;
transa[0] = 'N';
transb[0] = 'N';
m = ma->row; n = mb->column; k = ma->column;
alpha = 1.0; beta = 0.0;
a = (double*)malloc(sizeof(double)*m*k);
b = (double*)malloc(sizeof(double)*k*n);
c = (double*)malloc(sizeof(double)*m*n);
lda = m; ldb = k; ldc = m;
matrix2blas(ma,a,m,k);
matrix2blas(mb,b,k,n);
dgemm_(transa,transb,&m,&n,&k,&alpha,a,&lda,b,&ldb,&beta,c,&ldc);
blas2matrix(c,m,n,mc);
free(a);
free(b);
free(c);
return true;
}

BLAS の関数を C から呼び出すには、 blaswrap.h をインクルードする必要があるようです。従って、インクルードするファイルや、リンクするライブラリが書いてある、 vcblas.h を作成しました。内容は次の通りです。

#include <blaswrap.h>
#pragma comment(lib,"blas.lib")
#pragma comment(lib,"libF77.lib")
#pragma comment(lib,"libI77.lib")
#pragma comment(lib,"clapack.lib")

作成したプログラムが正確に動くかは、次のようにして確かめられます。

void c_dgemm_test() {
struct matrix a,b,c;
int m,n,k;
m = 3; n = 3; k = 3;
alloc_matrix(&a,m,k);
alloc_matrix(&b,k,n);
alloc_matrix(&c,m,n);
init_random();
petit_random_matrix(a,10,10); //ランダムな値代入
petit_random_matrix(b,10,10); //ランダムな値代入
show_matrix(a); //初期値表示
show_matrix(b); //初期値表示
multi_matrix_dgemm(&a,&b,&c);
show_matrix(c); //演算結果表示
multi_matrix(&c,a,b); //自前の行列積演算
show_matrix(c); //演算結果比較
free_matrix(&a);
free_matrix(&b);
free_matrix(&c);
}

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

A=
-1.00000  -4.00000  1.00000
5.00000  -4.00000  0.00000
5.00000  -2.00000  -1.00000
B=

0.00000  5.00000  5.00000
2.00000  -1.00000  1.00000
5.00000  -4.00000  -3.00000
Call dgemm_...
C=
-3.00000  -5.00000  -12.00000
-8.00000  29.00000  21.00000
-9.00000  31.00000  26.00000
Matrix multi...
C=
-3.00000  -5.00000  -12.00000
-8.00000  29.00000  21.00000
-9.00000  31.00000  26.00000

どうやら、 BLAS の呼び出しは成功しているらしいです。 BLAS のインストールも終了しました。では、 BLAS を使うと、どのくらい高速化されるか、試してみます。使う行列はもちろん、 1000 ×

1000 です。 BLAS の関数のプログラムは、先ほどのものを使いますので、最初に dgemm 以外の処理時間を計測しておきます。 10 回平均で、 381[ms] でした。あまり処理時間に影響を与えていないことが分かります。では、前々回作成したプログラムと実行速度を比較してみます。 CPU は Athlon900MHz です。

アルゴリズム使用時間 (ms)
自前・最適化なし  100922
自前・最適化あり  5810
標準 BLAS  34191

やはり、標準 BLAS では、マシンに特化されたプログラムにはかなわないようです。したがって、次に最適化 BLAS を導入してみます。

2.4 最適化 BLAS の導入

Athlon 用の最適化 BLAS は見つからなかったため、 ATLAS で作成されたとみられる、 Athlon 用の BLAS を使うことにしました。これは、 Netlib の ATLAS ページの Windows 用^{4}にあった、 atlas v3p3 WinNT ATHLON256L2.tgz を利用しました。これには、すでにコンパイル済みのものが入っているようです。では、 C から使ってみます。 ATLAS で作られた BLAS を参照するには、 cblaswr.c もインクルードする必要があるらしいです。また、 blas.lib の代わりに、 libatlas.a にリンクします。しかし、実際にリンクさせてみたところ、「libatlas.a(cblas xerbla.o) : error LNK2001: 外部 ” impure ptr” は未解決です」というエラーがでました。このエラーをなくそうと、いろいろ試しましたが、エラーをなくすことができませんでした。そこで代わりに、 Pentium �用の BLAS を使ってみたところ、すんなり動作しました。どうやら、 Athlon 用は ATLAS version3.0.3 で作られているのに対して、 Pentium �用は、 ATLAS version3.2.1 で作られていることが動いた原因ではないかと思います。ならば、 ATALS version3.2.1 を使って、自分でコンパイルすればよさそうです。しかし、今回は時間が無いので自分では作らず、 Pentium �用の BLAS で我慢します。また、 CBLAS を定義することで、標準 BLAS と最適化 BLAS とを変更できるようにするため、 vcblas.h を次のようにしました。

4http://www.netlib.no/netlib/atlas/archives/windows/

#include <blaswrap.h>
#ifndef CBLAS
#pragma comment(lib,"blas.lib")
#else
#include <cblaswr.c>
#pragma comment(lib,"libatlas.a")
#pragma comment(lib,"liblapack.a")
#pragma comment(lib,"libcblas.a")
#endif
#pragma comment(lib,"libF77.lib")
#pragma comment(lib,"libI77.lib")
#pragma comment(lib,"clapack.lib")

では、先ほどの 1000 × 1000 行列の積の計算と、所要時間を比べてみます。

アルゴリズム使用時間 (ms)
自前・最適化なし  100922
自前・最適化あり  5810
標準 BLAS  34191
最適化 BLAS(Pentium �用)  3444
最適化 BLAS(Pentium �用)  3395

自分でプログラム作るよりは、かなり速いです。もし Athlon 用の BLAS があれば、もう少し速くなるでしょう。数値計算を行うのなら、 BLAS を使ったほうが賢明な気がしました。

ここまでで、 CLAPACK と BLAS を C プログラムから使用する環境が整いました。次は、今回の課題を解いていきます。

3 問題の考察

3.1 問題 1

dgemm という関数を作るということでした。この関数は、次のような定義になっていました。

int dgemm_(char *transa, char *transb, integer *m, integer *
 n, integer *k, doublereal *alpha, doublereal *a, integer *lda,
 doublereal *b, integer *ldb, doublereal *beta, doublereal *c__,
 integer *ldc)

この関数は、行列積演算を行います。それは次のようなものです。

C = α × op(A) × op(B) + β × C

ここで、 α, β はスカラーで、 A, B, C は行列です。 op(A) は m × k 行列、 op(B) は k × n 行列、 C は m × n 行列です。また、個々の引数の細かい意味は、前のページを参考にしてください。

まずは、今回使うであろう、自分で定義した関数を確認します。

• void alloc matrix(struct matrix* m,int row,int col); row × column 行列の構築

• void free matrix(struct matrix* m); 行列の解放

• bool zero matrix(struct matrix m); 行列を 0 で初期化

• bool plus matrix(struct matrix* a,struct matrix m1,struct matrix m2); 行列 m1 と行列 m2 を足した結果を行列 a に格納

• bool times matrix(struct matrix* a,struct matrix m,double x); 行列 m をスカラー x で掛けた結果を、行列 a に格納

• bool transpose matrix(struct matrix* a,struct matrix m); 行列 m の転置行列を行列 a に格納

• bool copy matrix(struct matrix* a,struct matrix m); 行列 m の内容を行列 a にコピー

• bool blas2matrix(double* a,int m, int n, struct matrix* mt); BLAS での行列を、構造体の行列に変換

• bool matrix2blas(struct matrix* mt, double* a,int m, int n); 構造体の行列を BLAS での行列に変換

これらを使用して、プログラムを作成すると次のようになります。

int dgemm_(char *transa, char *transb, integer *m, integer *n,
        integer *k, doublereal *alpha, doublereal *a, integer *lda,
        doublereal *b, integer *ldb, doublereal *beta, doublereal *c__,

       integer *ldc)
{
struct matrix ma,mb,mc,mt,mt2;

//エラー

if ((*m < 0) || (*n < 0) || (*k < 0)) return -1;

//変数の初期化

alloc_matrix(&ma,-1,-1);alloc_matrix(&mb,-1,-1);
alloc_matrix(&mc,-1,-1);alloc_matrix(&mt,-1,-1);

//最初の次元？を設定

ma.row = *lda; mb.row = *ldb; mc.row = *ldc;
//op(A) の実行
alloc_matrix(&ma,*m,*k);
blas2matrix(a,ma.row,ma.column,&ma);
switch (transa[0]) {
case 'N': case 'n':

//そのまま

    break;
case 'T': case 't': case 'C': case 'c':

//転置行列

    transpose_matrix(&mt,ma);
    copy_matrix(&ma,mt);
    break;
default:
    free_matrix(&ma);free_matrix(&mb);free_matrix(&mc);free_matrix(&mt);
    return -1;
}
//op(B) の実行
alloc_matrix(&mb,*k,*n);
blas2matrix(b,mb.row,mb.column,&mb);
switch (transb[0]) {
case 'N': case 'n':

//そのまま

    break;
case 'T': case 't': case 'C': case 'c':

//転置行列

    transpose_matrix(&mt,mb);
    copy_matrix(&mb,mt);
    break;
default:
    free_matrix(&ma);free_matrix(&mb);free_matrix(&mc);free_matrix(&mt);
    return -1;
}

//行列積計算

if (ma.column != mb.row) {

//エラー

    free_matrix(&ma);free_matrix(&mb);free_matrix(&mc);free_matrix(&mt);
    return -1;
}
quick_multi_matrix(&mc,ma,mb,84);
//alpha 倍
if (!times_matrix(&mt,mc,*alpha)) {

//エラー

    free_matrix(&ma);free_matrix(&mb);free_matrix(&mc);free_matrix(&mt);
    return -1;
}
//C を beta 倍
if (*beta == 0) {
    copy_matrix(&mc,mt);
} else {
    free_matrix(&mc);alloc_matrix(&mc,mt.row,mt.column);
    blas2matrix(c__,mt.row,mt.column,&mc);
    alloc_matrix(&mt2,-1,-1);
    times_matrix(&mt2,mc,*beta);
    plus_matrix(&mc,mt,mt2);
    free_matrix(&mt2);
}
matrix2blas(&mc,c__,mc.row,mc.column);
free_matrix(&ma);free_matrix(&mb);free_matrix(&mc);free_matrix(&mt);
return 0;
}

まず最初に、与えられた引数が計算できるものであるか、確認します。次に、値の初期化を行います。ただ、 lda,ldb,ldc については、意味がよく分からなかったので、それぞれの行列の初期の値に使いました。ただ、意味はなしてません。その次に、 op(A) と op(B) を実行しています。とりあえず、行列の構造体に値を格納してから、転置行列にする場合のみ転置行列に変換しています。次に、前々回作成した、行列積を計算するルーチンを用いて、行列積を計算させます。その後、 alpha 倍したり、 beta が 0 でなければ C を beta 倍したりしています。最後にその結果を、元の配列に書き込んでいます。ただ、これだとオーバーベッドが大きくなりそうです。

まず、これが正常に動作しているのか、前で使った、 dgemm を呼び出すプログラムを利用して、今回作成したプログラムと、普通の行列積計算との結果の比較をしてみました。実行結果は次のようになりました。

A=
1.00000  5.00000  1.00000  -1.00000
3.00000  -3.00000  1.00000  3.00000
2.00000  4.00000  -2.00000  -4.00000

B=
1.00000  3.00000  3.00000
-3.00000  -3.00000  3.00000
-4.00000  -4.00000  0.00000
4.00000  -3.00000  5.00000
Call dgemm_...
C=
-22.00000  -13.00000  13.00000
20.00000  5.00000  15.00000
-18.00000  14.00000  -2.00000
Matrix multi...
C=
-22.00000  -13.00000  13.00000
20.00000  5.00000  15.00000
-18.00000  14.00000  -2.00000

どうやら、転置をしなかったり、定数倍しないときは正常に動作しています。また、 alpha の値を設定したり、 beta の値を 0 以外にして、 C の値に答えを加えるようにしたときも正常に動作しました。

次に、転置の設定がうまく動くか調べます。 transa と transb にそれぞれ’T’ を設定しました。また値を比較するため、自前の行列積関数に渡す行列は転置にしました。その結果、うまく動くことが確認できました。

A=
3.00000  -3.00000  0.00000  -1.00000
-1.00000  2.00000  -3.00000  4.00000
4.00000  -3.00000  0.00000  2.00000
B=
2.00000  -2.00000  1.00000
-3.00000  -4.00000  0.00000
-2.00000  0.00000  -3.00000
0.00000  -4.00000  0.00000
Call dgemm_...
C=
12.00000  -5.00000  -18.00000  4.00000
-13.00000  1.00000  15.00000  -8.00000
6.00000  12.00000  0.00000  12.00000
-8.00000  -13.00000  -4.00000  -16.00000
Matrix multi...

Transpose...
C=
12.00000  -5.00000  -18.00000  4.00000
-13.00000  1.00000  15.00000  -8.00000
6.00000  12.00000  0.00000  12.00000
-8.00000  -13.00000  -4.00000  -16.00000

最後に、 1000 × 1000 の行列を計算するのに、どれくらいかかるか計測しました。その結果、下のようになりました。

アルゴリズム使用時間 (ms)
行列積のみ  6142
自前 BLAS  7421

やはり、仮引数の形に合ったプログラムを書いたほうが。効率がよさそうです。

3.2 問題 2

Ax = b を解くための LAPACK の関数 dgesv (...) を作ります。今回は BLAS を利用するため、 BLAS に定義されている関数を確認しておきます。 dgesv は、

( \begin{matrix}a_{11} & a_{12} & ⋯ & a_{1n} \\ a_{21} & a_{22} & ⋯ & a_{1n} \\ ⋮ & ⋮ & ... & ⋮ \\ a_{n1} & a_{n2} & ⋯ & a_{nn} \end{matrix} ) ( x_{1} ; x_{2} ; ⋮ ; x_{n} ) = ( b_{1} ; b_{2} ; ⋮ ; b_{n} )

L = ( \begin{matrix}1 & 0 & ⋯ & 0 \\ l_{21} & 1 & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ l_{n1} & l_{n2} & ⋯ & 1 \end{matrix} ) U = ( \begin{matrix}u_{11} & u_{12} & ⋯ & u_{1n} \\ 0 & u_{22} & ⋯ & u_{1n} \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & u_{nn} \end{matrix} )

を解く関数です。したがって、 LEVEL2 BLAS を使うことになるでしょう。これは、連立一次方程式を解くので、三角行列に関係する BLAS を使うことになりそうです。なぜなら、 LU 分解を使うからです。 LU 分解とは、 A を下三角行列 L と上三角行列に分解して計算する方法です。

Ax = L(Ux) = b

において、 Ux = y とし、 Ly = b つまり、 y = L^{−1}b を解きます。次に、 Ux = y から、 x = U^{−1}y を求めることで解を得ます。この操作に一番あう BLAS の関数は、 dtrsv です。この内容は、 x ← A^{−1}x , x ← A^{−T} x (A : 三角行列) です。まさに、 LU 分解でのサブルーチンとなっています。

ではまず、この dtrsv の引数の意味を確認します。関数の形は次のようになっていました。

int dtrsv_(char *uplo, char *trans, char *diag, integer *n,
doublereal *a, integer *lda, doublereal *x, integer *incx)

• dtrsv DTRSV は、方程式を解きます。

Ax = b or A^{t}x = b

ここで、 b と x は n 個の要素をもつベクトルで、 A は n × n 行列であり、上か下の三角行列です。このルーチンには、きわめてまれにしか発生しないようなものはテストされていません。そのようなテストはこのルーチンを呼び出す前に実行されなければいけません。

• UPLO char[1] 型 UPLO は、マトリックスが上三角行列であるかした三角行列であるかを次に従って指示します。 UPLO = ’U’ or ’u’ A は上三角行列 UPLO = ’L’ or ’l’ A は下三角行列この引数は関数中で変化しません。

• TRANS char[1] 型 TRANS は、方程式を次のように解くことを指定します。 TRANS = ’N’ or ’n’, A*x = b TRANS = ’T’ or ’t’, A’*x = b TRANS = ’C’ or ’c’, A’*x = b この引数は関数中で変化しません。 • DIAG char[1] 型

DIAG は、 A が三角行列であるか次に従って指定します。 DIAG = ’U’ or ’u’ A は三角行列であると仮定 DIAG = ’N’ or ’n’ A は三角行列であると仮定しないこの引数は関数中で変化しません。

• N int 型 N は、行列 A の大きさを指定します。 N は少なくとも 0 でなくてはなりません。この引数は関数中で変化しません。

• A double 型配列 (LDA,n) 次元 UPLO = ’U’ or ’u’ を指定する前に、上三角行列の配列 A の (n,n) が上三角行列に含まれている必要があり、下三角の部分については厳密に言及しません。 UPLO = ’L’ or ’l’ を指定する前に、下三角行列の配列 A の (n,n) が下三角行列に含まれている必要があり、上三角の部分については厳密に言及しない。 DIAG = ’U’ or ’u’ のとき、 A の対角線の要素は、言及されません。ただし、 1 であることを仮定されています。この引数は関数中で変化しません。

• LDA int 型 LDA は、（サブ）プログラムを呼ぶときに宣言した A の最初の次元を指定します。 LDA は少なくとも max(1,n) でなければなりません。この引数はは関数中で変化しません。

• X double 型配列 (1+(n-1)*abs(INCX)) 次元以上増加された配列 X は、ベクトル b の左側に n の要素を含まなければなりません。終了時に X はベクトル x の結果として上書きされます。

• INCX int 型 INCX は、 X の要素の増加分を指定します。 0 であってはなりません。この引数は関数中で変化しません。

• LEVEL2 BLAS の関数です。

次に、 LAPACK の関数、 dgesv の関数定義について確認します。関数の形は次のようになっています。

int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
*lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)

• dgesv DGESV は、線型方程式のリアルシステムでの解答を計算します。

A × X = B

ここで、 A は n × n 行列で、 X と B は、 n × nrhs 行列です。部分ピボットと行交換をを含む LU 分解は、 A を次のように分解することで得られます。

A = P × L × U

ここで、 P は行列の交換であり、 L は下三角行列、 U は上三角行列です。 A を形式的に分解することは、 A × X = B の方程式を解くのに使います。

• N int 型線形方程式の値です。行列 A の大きさを指定します。 N は 0 以上でなくてはなりません。

• NHRS int 型左端の値です。行列 B の列数です。 NHRS は 0 より大きくなければなりません。

• A double 型配列 (LDA,N) 次元行列 A の大きさは、 n × n です。終了時、 A = P × L × U の分解の仕方で、 L と U に分解されています。 L の対角線の要素は保存されていません。

• LDA int 型行列 A の次元から導き出されます。 LDA は max(1,N) 以上でなければなりません。

• IPIV int 型配列 n 次元ピボット選択による行列 P の交換の定義をします。行列の i 行目は IPIV(i) 行目だったものです。

• B double 型配列 (LDB,NRHS) 次元行列 B の左端は (N × NHRS) です。終了時に、 info = 0 ならば、解を解くのに成功しています。

• LDB int 型配列 B の最初の次元です。 LDB は max(1,N) 以上でなければなりません。

• INFO int 型 = 0 演算は成功 < 0 INFO = -i とすると、 i 番目の引数が、不正な値 > 0 INFO = i とすると、 U(i,i) がちょうど 0 である。分解は成功している。しかし、 U はまれに見る行列であり、答えを計算できなかった。

では、プログラムを作成します。今回作成するものは、基本的に LU 分解であるため、プログラムの細かい説明はしません。前回のレポートの内容を、少し手直ししただけです。プログラムは、次のようになりました。

int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
*lda, integer *ipiv, doublereal *b, integer *ldb, integer *info)
{
int i,j,k,p;
int N;
double x,r;
/* BLAS の引数 */
char uplo[1]; char trans[1]; char diag[1]; double* b_i; int incx;
N = *n;

/***** エラー *****/
if (N < 0) {*info = -1;return -1;}
if (*nrhs < 0) {*info = -2;return -1;}
if (*lda < 1 || *lda < N) {*info = -4;return -1;}
if (*ldb < 1 || *ldb < N) {*info = -7;return -1;}
/***** LU 分解 *****/
for (k=0;k<N;k++) ipiv[k] = k; /* ピボット交換のベクトルを初期化 */
for (k=0;k<N-1;k++) {
    for (i=k+1;i<N;i++) { /* 第 k 列より下を 0 にする */
        if (a[k+k*N] == 0) { /* 部分ピボット選択 */
            p = k; x = dabs(a[k+k*N]);
            for (j=k+1;j<N;j++) { /* 絶対値が最大の行を検索 */
                if (x < dabs(a[j+k*N])) {
                    p = j; x = dabs(a[j+k*N]);
                }
            }
            if (x == 0) {*info = -k;return -1;}
            for (j=k;j<N;j++) { /* k 行と p 行の入れ替え */
                x = a[k+j*N]; a[k+j*N] = a[p+j*N]; a[p+j*N] = x;
            }
            ipiv[k] = p;
        }
        r = a[i+k*N] / a[k+k*N];
        a[i+k*N] = r;
        if (r == 0) continue; /* 計算の必要なし */
        for (j=k+1;j<N;j++) { /* 行列部分の演算 */
            a[i+j*N] -= r*(a[k+j*N]);
        }
    }
}
for (i=0;i<N;i++) { /* LU 分解の結果が計算できるものか */
    if (a[i+i*N] == 0) {
        *info = -i; return -1;
    }
}
/***** LU 分解完了 *****/

/***** LU 分解の結果を用いて連立一次方程式をとく *****/

b_i = (double*)malloc(sizeof(double)*N);
uplo[0] = 'U'; trans[0] = 'N'; diag[0] = 'U'; incx = 1;

    for (k=0;k<*nrhs;k++) {
        for (i=0;i<N;i++) b_i[i] = b[i+k*N]; /* b から一列抜き出す */
        for (i=0;i<N;i++) { /* ピボット選択 */
            if (i != ipiv[i]) {
                x = b_i[i]; b_i[i] = b_i[ipiv[i]]; b_i[ipiv[i]] = x;
            }
        }
        /* L を解く */
        uplo[0] = 'L'; diag[0] = 'U'
        dtrsv_(uplo,trans,diag,n,a,n,b_i,&incx);
        /* U を解く */
        uplo[0] = 'U'; diag[0] = 'N';
        dtrsv_(uplo,trans,diag,n,a,n,b_i,&incx);
        for (i=0;i<N;i++) { /* ピボット選択 */
            if (i != ipiv[i]) {
                x = b_i[i]; b_i[i] = b_i[ipiv[i]]; b_i[ipiv[i]] = x;
            }
        }
        for (i=0;i<N;i++) b[i+k*N] = b_i[i]; /* b へ一列戻す */
    }
    free(b_i);
    *info = 0; return 0;
}

行列 A を、 LU に分解するところまでは、前回の LU 分解のプログラムと同じです。 LU を用いて、値を求めるときに、 BLAS の関数を利用しています。ここではまず、行列から答えを求めるベクトルを抜き出します。そして、ピボット選択の結果にあわせて、ベクトルの値を入れ替えています。次に、 Ly = b の値を、 BLAS で解きます。このとき L は下三角行列であり、対角線の要素は全て 1 であることから、 uplo = ’L’、 diag = ’U’ としてあります。そして、 dtrsv に計算させると、 b i には、 y = L^{−1}b の結果が格納されています。次に、同様にして Ux = y を計算します。 U は上三角行列で、対角線の要素はある値になっているため、 uplo = ’U’、 diag = ’N’ と指定します。最後に、 b i の並びをピボット選択する前に直し、値を b に代入します。この操作を、 nrhs で指定された分だけ行えば、この関数は完成です。

では、実際にうまく動くのか、前回作成した LU 分解のプログラムと比較してみます。

N = 3
Called dgesv_
4.500000
-3.000000
2.500000
4.50000

-3.00000
2.50000

また、 BLAS を入れ替えたら、本当に高速化されるのか確かめてみます。予想としては、全く高速化されないと思われます。なぜなら、最初の行列 A に依存する計算の回数は、 n^{3}/3 + O(n^{2}) であるのに対して、 BLAS の担当しているベクトル b に依存する計算の回数は、 O(n^{2}) だからです。したがって、 BLAS が高速化されてもさほど速度に影響がでそうもないということです。では、実際に試してみます。試した大きさは、 1000 個の連立一次方程式を解くものです。

アルゴリズム使用時間 (ms)
自前 dgesv + 標準 BLAS  48976
自前 dgesv + 最適化 BLAS  48970
CLAPACK dgesv + 標準 BLAS  10154
CLAPACK dgesv + 最適化 BLAS  1419

予想通りです。どうやら、 BLAS 高速化によってこの関数が高速化されるためには、前半の LU 分解を BLAS で書く必要があるようです。また、 CLAPACK の関数は最適 BLAS により高速化されていることから、 LU 分解のプログラムも、 BLAS で書かれていると思われます。

4 感想

やはり、 CLAPACK の関数は性能がいいです。さすが、プロが書いたものだけあります。また、 ATLAS で自動生成された BLAS の関数も、とても性能がいいです。自分で書いたプログラムより、数倍から数十倍速くなります。

この結果から分かることは、利用できるライブラリは積極的に利用したほうがいいということです。やはりこれからは、いかに多くの API やライブラリを知っていて、それを使いこなせるのかがより重要になります。

しかし、今回の課題も大変でした。さて、どれくらいこのレポートに時間を使ったでしょうか？少なくとも一桁ではありません。ただ単に、僕の理解する速度が遅いだけかもしれませんが。これからまた、このように大変なレポートが出題されないことを願うのみです。