Contents

  • ベクトルの格納
  • 行列の格納
  • 行列のパックド形式
  • 複素数の場合
  • 呼び出し方
  • コンパイル
  • incX,ldAについて
  • 単位三角行列の指定について
  • 注意点
  • CBLASを使いたい！
  • なんで、BLASを直接リンクする方針なのですか

ベクトルの格納

BLASにおけるベクトルの格納は、単純な配列になっています。例えば、(1.0,2.0,3.0)と値が入っているベクトルは

double x[3] = {1.0, 2.0, 3.0 };

のように書けます。以降これを列ベクトル（つまり縦に並んだ物）として扱っていくことになります。行ベクトルを扱いたい場合は、扱う計算の方のベクトルを全部転置してコーディングしてください。

行列とベクトルの積を扱いたいという場合も、ベクトルは右から係る形で計算されます。

行列の格納

このサイトでは、FORTRANのルーチンを直に呼び出すやり方を説明しているので、行列の格納形式には注意が必要です。不正確ですが、結果から申し上げると、与えた行列が転置されます。

double A[3][3]={
{1.0, 2.0, 3.0},
{4.0, 5.0, 6.0},
{7.0, 8.0, 9.0},
};

と、
1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 9.0

なる行列を与えたとします。このとき、この行列をルーチンに渡したとします。すると、
1.0 4.0 7.0
2.0 5.0 8.0
3.0 6.0 9.0

という行列を与えたことになっています。つまり、転置されたものが渡っているように見えるわけです。実際、dgemv_（一般行列とベクトルの積）に先程の、xとAを渡して、Axを計算させると、
30.0
36.0
42.0
が返ってきて、期待していた
14.0
32.0
50.0
とはなっていません。

個人的には、double A[3*3]などと、行列の全ての要素数を一次元配列で用意して(i,j)成分のアクセスには、A[i+j*3]のようにアクセスすることをおすすめします。 
#define MatA(i,j) (A[(i)+(j)*3])
のようなマクロを用意してやると良いと思います。MatA(1,2)のようにアクセス可能です。（当然ですが、要素は0〜2になりますので、この場合はAの2行3列にアクセスしていることになります）
ちなみに、長方行列 3行4列の場合ではdouble A[3*4]となって、A[i+j*3]となります。
このようにアクセスすると、上記の問題は起こりません。（もともと転置した場所をアクセスするようになるためです。というか、FORTRANではこういう並びになってるんです）また、動的にメモリを確保する際にも同様に取り扱うことができ、非常に便利です。

なお、このようなマクロの使用時は、引数の扱いに十分注意してください。「関数マクロ　危険」などで検索すると実例は多く出てくると思われます

行列のパックド形式

対称行列（エルミート行列）や三角行列など、行列の上半分のみや下半分のみが必要な行列については、メモリの使用量を圧縮する形式が用意されています。

この形式で行列を渡すことができるルーチンでは、引数uploで'U'か'L'を指定し、上三角部分（右上の部分）か下三角部分（左下の部分）のどちらを用いているかを指定します。
'U'を指定した場合、i<=jの部分をアクセスすることになりA[(i+j(j-1)/2)]にA(i,j)成分が格納され
'L'を指定した場合、i>=jの部分をアクセスすることになりA[(i+(2*n-j)*(j-1)/2)]にA(i,j)成分が格納され
 メモリの使用量を半分にすることが出来ます。(nは行列の行または列：今対称、三角行列を扱うので正方行列になっています）
たとえば、以下のような対称行列を格納することを考えると
1.0 2.0 3.0
2.0 5.0 6.0
3.0 6.0 9.0

パックド形式は上半分を格納する場合メモリには
1.0 2.0 5.0 3.0 6.0 9.0
と順番に格納されているばあい以上のように認識されます。

複素数の場合

複素数を扱う場合は、
struct { double r, i; }
のような構造体を用意する必要があります。f2c.hというFORTRANのルーチンをC言語から取り扱うためのヘッダファイルでは、これをdoublecomplex（単精度(float)版はcomplex）と定義しています。上記のdoubleの代わりにこれを用いて下さい。

例えば、(1.0+4.0i, 2.0+5.0i, 3.0+6.0i)と値が入っているベクトルは
doublecomplex z[3] = {{1.0, 4.0}, {2.0, 5.0}, {3.0, 6.0}};
と書くことが出来ます。(C言語では構造体の初期化も書けますので）
この場合、値へのアクセスは、
z[0].r
z[0].i
のようにして行うことが出来ます。
当然ですが、同様の構造体を自前で用意しても構いません。
blas.hでも、complex,doublecomplexのように定義しています。

呼び出し方

このサイトで紹介しているルーチン名をそのまま書けばOKです。が、特にC++で用いる場合には、プロトタイプ宣言が必要になります。CBLASにはヘッダファイルがありますので、これを読めばOKなのですが、このサイトではこれを用いませんので、プロトタイプ宣言を自分で書かなくてはいけません。(Cでは無くても通りますが)

それは面倒だという人のために、blas.hを作成しておきました。ここにはこのサイトで紹介しているルーチンの全てのプロトタイプ宣言が書かれています

たとえば、倍精度実数ベクトルのスカラ倍、dscalを用いる場合、プロトタイプ宣言は
int dscal_(int* N, double* alpha, double* X, int* incX);
となります。最後にアンダースコア_が付いていることに注意してください。C++で使う場合は、extern "C"で囲む必要があります。

とは言っても、それほど難しいことはありませんし、一度作ってしまえばさほど手間ではありません。

この宣言を見ていただければ分かると思いますが、全ての引数が参照渡し、ポインタで渡す事になっています。これはFORTRANの仕様によるものです。そのため、これを呼び出す際には、すべての値を一度変数につめ込みアドレスを渡す必要があります。例えば、alpha=2として、ベクトルの要素を2倍にしたい場合では
double alpha=1.0,;
としてから、&alphaなどとしなくてはなりません。ただし、文字を渡す際には文字列としてそのまま書いても構いません。Cにおいて、文字列リテラルは全てポインタになるからです。

例に挙げたdscal_は以下のように呼び出せます。

double x[3] = {1.0, 2.0, 3.0 };
int n=3;
double alpha=2.0;
int one=1;

dscal_(&n, &alpha, x, &one);

サンプルコードも参考にしてください。

コンパイル

ソースコード本体とBLAS、FORTRANのライブラリをリンクさせればOKです。コンパイルは環境によって異なりますので参考程度にしていただきたいのですが、Linux環境でGCCを用いてGotoBLAS2をリンクさせる場合
gcc source.c -lgoto2 -lgfortran
とすることになります。（パスの通っているところにlibgoto2.aが置かれている場合） パスが通っていない場合は
gcc source.c -L/home/username/lib/goto2 -lgoto2 -lgfortran
などのように、-Lでライブラリが置いてあるディレクトリを指定する必要があります
サンプルコードも参考にしてください。

incX,ldAについて

リファレンスには、incXはインクリメント幅などと書いてあり意味が分からなくなるのですが、incX,incYなどは1を、ldA,ldBなどは、行の数を指定するのが通常です。通常はこれらを引数に指定するようにしてください。詳しくは、サンプル(DGEMVの事例)で例を用いて説明しています。

単位三角行列について

三角行列を用いる場合、diagという引数が付いてきて、単位三角行列かどうかを指定することが出来ます。これは、"N"の時は、対角成分をそのまま読むが、"U"の時は、対角成分を無視して1として扱うという意味になります。なぜこのような指定があるかというと、LU分解などを行った際には、下三角行列部の対角成分は1になりますが、これは格納されません。このような場合にdiag="U"を指定することになります。

注意点

BLASのほとんどのルーチンは、与えた引数の行列やベクトルの少なくともひとつが破壊されることに注意してください。計算前の値を用いる必要がある場合は、与えた引数が破壊されていないかどうか必ず確かめるようにしてください。破壊される場合は必要に応じて?copy_で予めコピーを取っておきましょう。

また、BLASには疎行列を扱うルーチンはありません　SparseBLASといったルーチンがありますのでそちらを使用することになります。

CBLASを使いたい！

CBLASはC言語向けのBLASのインタフェースです。これを用いたいという場合は、CBLASをBLASとあわせてインストールし、cblas.hをインクルードしてください。　基本的な関数の使い方に違いはありませんが（内部で同じものを呼んでいるので）、以下のような違いがあります。
  • 名称が例えば、cblas_dgemmなどのようになる
  • 引数に値渡しが許される
  • 行列の列優先、行優先などの設定ができる
  • 転置などの設定が文字列ではなく列挙型
これらは、cblas.hのプロトタイプ宣言を見れば、BLASの利用法さえ理解していれば十分に利用が可能です。

なんで、BLASを直接リンクする方針なのですか

筆者の趣味です。その大きな理由は、LAPACKにあります。LAPACKをCで用いるためのライブラリとして、CLAPACKという物がありますが、これはインタフェースではなく、f2cというソフトウェアで機械的にC言語に翻訳された（というなんとも乱暴な）代物です。正直、信頼性がかなりイマイチであると考えています。で、どうせLAPACKはFORTRANのものを直接呼ぶことにするなら、BLASもFORTRANのものでいいじゃないか。という考えに至りました。以上で見るように、直接FORTRANのライブラリをリンクする場合には、行列の格納などについて注意をはらう必要がありますが、実は、CBLASではそこまで注意が必要ではありません。それなら、CBLASを使いたいのですが、結局のところ、LAPACKを使う以上はその注意をしなくてはいけないわけで、CBLASを使うメリットがほぼ無くなります。また、両方を混在させることに依る混乱も避けたいところです。ということで、全てFORTRANを直接呼ぶことにしました。CBLASで毎回指定しなくてはいけないCblasRowMajorなども書かなくて済むので、ソースもすっきりしましたので、そこまで悪いことはないと思います。逆に言えば、BLASだけしか、使わないという方はCBLASを使うほうが楽になる点はあるかと思います。そのあたりは自身で判断してください。