ATLASを使う - Risky Dune

行列演算がボトルネックのコードを書いたので, Automatically Tuned Linear Algebra Software (ATLAS)を導入した.

キャッシュサイズなどを考慮して行列演算を行なってくれるため, 自分で適当に書いた奴より速い. uBlasより速いか検討してないが, まあ多分速いだろう. ただ, インターフェイスとしてはおそらくuBlasの方が優れているのでそちらの使用をまず検討した方がいいかもしれない.

インストール

環境はUbuntu12.04です.

手順を列挙すると

aptからgfortranを入れる
CPU frequency scalingを切る
ATLASをコンパイル & インストール

おそらくgfortranはATLASがFortran上に構築されてるから必要.
CPU frequency scalingは計算量があまり必要ないときに動作周波数を抑えて, 消費電力を小さくする技術.
最適化の上でこれが邪魔になるらしい.
Ubuntu上でこれを切るにはrootで

echo performance > /sys/devices/system/cpu/cpu0/cpufreq/scaling_governor
echo performance > /sys/devices/system/cpu/cpu1/cpufreq/scaling_governor

みたいに論理cpuごとにscaling_governor(たぶん周波数調整の戦略)を書き換える.

あとはインストール. だいたいATLAS公式のインストールガイドとATLASに付属するINSTALL.txtを見ればわかる. 以下はこちらの環境での手順.

まずatlas本体とlapackをダウンロード.
今回使ったのはlapack-3.4.2.tgzとatlas3.10.0.tar.bz2. 以下, 記憶に従って書いてるので間違ってたらごめん.

cd ATLAS
mkdir my_build ; cd my_bulid
../configure --shared --with-netlib-lapack-tarfile=<path to lapack tarfile>
make

sharedは共有ライブラリが欲しい人用. 静的ライブラリだけでいいならいらない. --with-netlib-lapack-tarfileはlapackのインターフェイスを提供して欲しい人用. これが無いと線形方程式を解くインターフェイスとかは提供してもらえない. lapackのファイルはtgzのままで確かいけたと思う.

makeした後の配置はmake installで多分できる. my_build/libに出来たライブラリとATLAS/includeのincludeファイルをコピーしてもいいけど. ただ, コンパイル過程でのみ必要なファイルも混ざってるので注意.

C++で使う前に

ATLASはCとしてコンパイルされているのでC++から呼ぼうとするとundefinedうんたらかんたらと言われてしまいます.
そのためヘッダファイルに少し手を入れる必要があります.

#ifdef __cplusplus
extern "C" {
#endif

と対応する閉じ括弧をcblas.hとclapack.hに前もって追加してください.

使い方 (ちょっと怪しい)

BLASはAPIの名前なので, 複数のBLAS実装が存在します(GotoBLASとかATLASとか). そして使い方は同じであるため, BLASの使い方について調べればそれはATLASにも適応できる筈！今回はCBLASというC用のインターフェイスを用います.
CBLASの真っ当な日本語ドキュメントがちょっと見当たらない. こちらのページはかなり良くまとまっていますが, 直接Fortranとやりとりするのでちょっと面倒.
関数名までわかっているならMacのリファレンスがわりと親切.

とりあえず簡単な使い方と, BLASの名前付け規約について説明します.
行列とベクトルの積を行うcblas_?gemvを例に見て行きます. ?の部分はs,d,c,zのどれかが入り

s : 単精度
d : 倍精度
c : 単精度複素数
z : 倍精度複素数

を対象にします. 複素数関係は使った事無いです. 今回はcblas_sgemvを使います.

cblas_sgemvは正確には次の計算を行います

$y = \alpha * Ax + \beta * y$

計算に使ったyに結果を代入するのに注意してください. 「俺はただ行列積をしたいだけなのに, 余計なalphaとかyを導入しやがって」と思うかもしれません. 僕も思います. ただ, ちゃんとソースを読んでいないのでわかりませんが, alpha=1.0のときやbeta=0.0のときに条件分岐があるっぽいので, そこで余計な演算をしないようになっている「かも」しれないです.

次がsgemvの定義です.

void cblas_sgemv (
const enum CBLAS_ORDER Order,
const enum CBLAS_TRANSPOSE TransA,
const int M,
const int N,
const float alpha,
const float *A,
const int lda,
const float *X,
const int incX,
const float beta,
float *Y,
const int incY
);

それぞれの意味は次のようになります. 後に載せるサンプルと見比べてください.

CBLAS_ORDER : AがRow MajorかCol Majorか. Row MajorならCblasRowMajorをセット
TransA : Aを転置(と共役)して用いるかどうか. とくにいらないならCblasNoTrans
M : Aの行 (転置したらNと入れ替わる)
N : Aの列 (転置したらMと入れ替わる)
lda : Row Majorなら列数. Col Majorなら行数.
incX : Xのいくつおきの要素を使うか. 普通は1
incY : 同上

これを使ったサンプルが以下のとおりです. Row Major Orderの例になっています
blas_blog.cpp

#include <stdio.h>
#include <cblas.h>
#include <clapack.h>

int main(){
  /*
    A = 1 2 3
        4 5 6
   */
  float A[6] = {1,2,3,4,5,6};

  /*
    x = 1
        2
        3

    y = 4
        5
   */
  float x[3] = {1,2,3};
  float y[2] = {4,5};

  cblas_sgemv(CblasRowMajor,
              CblasNoTrans,
              2,  /*M*/
              3,  /*N*/
              1.0, /*alpha*/
              A,
              3, /*lda*/
              x,
              1, /*incX*/
              1.0, /*beta*/
              y,
              1);

  for(int i=0; i<2; i++){
    printf("%.1f ", y[i]);
  }
  printf("\n");
}

clapack.hは今回不要ですが, lapack系のメソッドを使うときはいります.

Makefile

CC=g++ 
OBJS=blas_blog.o
INCS=-I/usr/local/atlas/include
LIBS=-L/usr/local/atlas/lib -llapack -lcblas -latlas

.cpp.o:
	$(CC) -c $(INCS) $< 

main:$(OBJS)
	$(CC) -o $@ $<   $(LIBS)

ライブラリの指定順序は意味を持つので注意してください. 実行結果は次のようになるはずです.

18.0 37.0

numpy/scipyと一緒に使う場合

numpy/scipyのバックエンドとしてATLASを使うことができます.
やり方は忘れましたが, site.cfgに

[DEFAULT]
library_dirs = /usr/local/lib:/usr/local/atlas/lib
include_dirs = /usr/local/include

と設定してやらないとATLASを見つけてくれなかったりしました.

本題なのですが, このバックエンドにATLASを使った場合, cpu frequency scalingをONにしているとデフォルトの行列演算よりも大分遅くなる現象が確認されました. 使っているときはcpu frequency scalingを切った方がいいかもしれません.