kashiの日記

検索条件

全5件 (1/1ページ) 記事本文を表示

2015/06/29(月)kv-0.4.22

kvライブラリを0.4.22にアップデートしました。

今回は、defint_power_r, defint_log_r, defint_power_autostep_r, defint_log_autostep_rの4つの端点特異性を持つ積分の関数のbug fixです。これら右側の端点特異性を扱う関数は、時間を反転させて左側versionに任せているのですが、その時間反転部分にバグがありました。これらの関数を使う例題を作っていなかったので、気付きませんでした。

その他、細かいところをいくつか修正しています。

コメント（0件）

2015/06/28(日)octave+openblas+intlabで丸め変更を検証

精度保証

INTLAB version 9をubuntu上のoctaveで使うの続きです。このopenblasはパッケージで入れたものなので、丸めの向きの変更が正しく効くかどうか分かりません。そこで、行列積について、丸めの向きの変更が効いているかどうかテストしてみました。前の記事で入れたubuntuと、ついでにwindowsでも試してみました。

行列積はn³回の乗算とn²(n-1)回の加算を含んでいますが、これらの全ての演算について、上向き丸めと下向き丸めで結果が異なる筈であるような入力値を与え、計算結果が同一になる(すなわち丸め変更が効いていない)ことが発生することがないかどうか試すプログラムを作ってみました。matlabのプログラムは書いたことがないのでこれが普通の書き方かどうかはよく分かりません。

function testmm(n)
	disp('testing multiplication...')
	flag = 1;
	for k=1:n
		a=zeros(n);
		a(:,k)=ones(n,1)*1/3;
		b=zeros(n);
		b(k,:)=ones(1,n)*1/3;
		c = (rad(intval(a) * b) ~= 0) + 0.0;
		if prod(prod(c)) == 0
			disp('multiplication error!')
			flag = 0;
			break
		end
	end
	if flag == 1
		disp('multiplication OK!')
	end

	disp('testing addition...')
	flag = 1;
	for k=2:n
		a=zeros(n);
		a(:,1)=ones(n,1);
		a(:,k)=ones(n,1)*2^(-27);
		b=zeros(n);
		b(1,:)=ones(1,n);
		b(k,:)=ones(1,n)*2^(-27);
		c = (rad(intval(a) * b) ~= 0) + 0.0;
		if prod(prod(c)) == 0
			disp('addition error!')
			flag = 0;
			break
		end
	end
	if flag == 1
		disp('addition OK!')
	end
end

これをtestmm.mという名前で作成し、testmm(100)などとintlab上で実行してみました。なお、数値は行列サイズで、このプログラムはn⁴に比例した時間がかかるのであまり大きな行列ではテストできません。実行結果は、次のようになりました。

	ubuntu 14.04	windows + octave 3.8.2	windows + octave 4.0.0
シングルコア	OK	OK	OK
マルチコア	NG	OK	NG

実行環境は次の通り。

ubuntu 14.04	前記事で書いたubuntu 14.04 + パッケージで入れたoctave, openblas
windows + octave 3.8.2	windows7 + INTLAB公式のoctaveページからダウンロードしたoctave 3.8.2
windows + octave 4.0.0	windows7 + Octave公式ページからダウンロードしたoctave 4.0.0

さて、この結果はある程度予想できるものでした。これら3つの環境では全て内部でopenblasが用いられています。openblasのソースファイルのMakefile.ruleというファイルの中に、

# CONSISTENT_FPCSR = 1

という行があります。これは、マルチスレッド環境において浮動小数点数のControl/Statusレジスタをコピーするか否かを決めていて、ここをuncommentするとopenblasが起動されたときの丸めの向きの状態がマルチスレッド時の全てのスレッドに伝わるようになります。この結果から、ubuntuのパッケージとoctave公式のopenblasはここがコメントのままmakeされていて、INTLABページのものはここのコメントが外されてmakeされていると推定できます。

ここからは、ubuntu 14.04の場合について、openblasのパッケージを作りなおしてマルチコア環境でも丸めの向きの変更が全てのスレッドに伝わるようにします。野良ビルドでも良いのですが、ubuntuのパッケージシステムとなるべく合わせた形で作ってみました。

sudo apt-get install devscripts
適当な作業directoryを作ってそこに移動
sudo apt-get source openblas でソースコードを取得。
sudo apt-get build-dep openblas でmakeするのに必要なものをインストール。
cd openblas-0.2.8/ でソースdirectoryに入る。
Makefile.ruleを編集し、コメントアウトされている
```
CONSISTENT_FPCSR = 1
```
を有効化する。
dch -i で、version番号を例えば0.2.8-6ubuntu2のように進め、changelogを適当に書く。
dpkg-buildpackage -b -uc -us でmake。
これで必要なdebファイルが生成されたので、cd .. で一つ上に上がり、
sudo dpkg -i libopenblas-base_0.2.8-6ubuntu2_amd64.deb
sudo dpkg -i libopenblas-dev_0.2.8-6ubuntu2_amd64.deb

別のマシンで入れ替え作業を行う場合は、最後に使った2つのdebファイルだけあればOKです。一応作成したdebファイルを上げておきますが、自分で作成することをお勧めします。

libopenblas-ubuntu-rounding.zip

ところで、INTLAB公式のoctaveページからダウンロードしたoctave-3.8.2-installer.exe、Aviraさんが

みたいに反応するんだが、大丈夫だろうか (自分は仮想環境で実行してすぐに消しました)。Aviraさんはときどき過剰に反応するからそのせいだと信じよう。

コメント（0件）

2015/06/24(水)INTLAB version 9をubuntu上のoctaveで使う

精度保証

INTLAB version 9の重要性の記事で書いたように、INTLAB version 9を使えばmatlabに頼らずにオープンソースな環境で精度保証付き数値計算を行うことができます。今回、ubuntu 14.04にoctaveをインストールし、その上でINTLABを動かしてみたので、その様子を記録しておきます。

INTLABは、http://www.ti3.tu-harburg.de/rump/intlab/で、50ユーロで購入できます。有償なのが残念ですが、PayPal経由で支払いを済ませると一時間だけ有効なダウンロードリンクがメールで送られてきます。ダウンロードしたファイルは完全にソースコードを見ることができます。また、Linux/Mac/Windows用、あるいはmatlab用/octave用といった区別はなく、すべてこのファイルでOKです。

octaveのインストールは、ubuntu 14.04で、

sudo apt-get install octave
sudo apt-get install octave-doc
sudo apt-get install octave-htmldoc
sudo apt-get install liboctave-dev

のように行いました。特に最後のliboctave-devが必要なことに注意してください。また、ubuntu 14.04でBLASを使うの記事のようにATLASやopenblasを導入すれば、octaveから使われるBLASがそれに置き換わります。openblasを入れることをお勧めします。

ダウンロードした INTLAB.zip を任意のdirectoryに展開します。例えば ~/Intlab_V9/ としましょう。まず、octaveを起動し、必要ならば cd Intlab_V9 などとしてINTLAB directoryに移動し、

startintlab

と実行します。これは、そこのdirectoryにある startintlab.m を実行しています。うまく行けば、これだけでintlabが使える状態になる筈です。初回起動時は、裏で

Intlab_V9/setround/octave/setround.mex
Intlab_V9/INTLAB_Data_Intlab_Version_9_x86_64-pc-linux-gnu3_8_1.mat

の2つのファイルが作られました。1つめは丸めの向きを変えるためのファイルで、2つめは数学関数の誤差を調査し記録しておくためのファイルです。1つめのファイルを作成するために、先に入れた liboctave-dev が必要になります。C++コンパイラもインストールされている必要があると思います。2つめのファイルは内部にdirectory情報を含んでいるらしく、intlab一式の場所を移してしまうと動かなくなります。その場合はそのファイルを削除して再度 startintlab を実行すれば大丈夫です。

実際に区間演算をしてみます。T先生からの挑戦状で計算した例題です。

octave:1> format long
octave:2> x=intval(1);
octave:3> y=intval(1);
octave:4> for i=2:10; z=(1+2*y)/(x*y^2); x=y; y=z; i, z, end
i =  2
intval z = 
   3.00000000000000
i =  3
intval z = 
   0.77777777777777
i =  4
intval z = 
   1.40816326530612
i =  5
intval z = 
   2.47448015122873
i =  6
intval z = 
   0.68995395922102
i =  7
intval z = 
   2.0203934747424_
i =  8
intval z = 
   1.7898074714308_
i =  9
intval z = 
   0.7075878600520_
i =  10
intval z = 
   2.69514211796___
octave:5>

のように計算出来ました。最初の intval を affari に変更すればaffine arithmeticで計算させることもできます。

また、行列積を計算させてみると、

octave:1> a=rand(5000);b=rand(5000);
octave:2> tic;a*b;toc
Elapsed time is 3.20024 seconds.
octave:3> tic;intval(a)*b;toc
Elapsed time is 8.66678 seconds.

くらいでした。core i7 4770Kですがvmwareを挟んでいるので生で使うと更に速いでしょう。

なお、octave起動時に毎回startintlabを実行させるのは面倒なので常時実行させるには、

/usr/share/octave/site/m/startup/octaverc
home directoryの.octaverc
current directoryの.octaverc

のいずれかのファイルに、

addpath('/intlabを入れたパス/Intlab_V9')
startintlab

と書いておけば良さそうです。

余談ですが、octave起動時に

octave --force-gui

とすると、

のような感じで起動します。変数一覧とか見られて便利そう。

openblasで丸め変更に問題はないかとか、検証する必要はあるかと思います。

自分はmatlab嫌いなものでINTLABの使い方はよく分かっていませんが、お楽しみ下さい。

コメント（0件）

2015/06/21(日)ubuntu 14.04でBLASを使う

ubuntu

ubuntu 14.04上でBLASをパッケージで入れて使ってみたのでメモを残しておきます。

BLASは、Basic Linear Algebra Subprogramsの略で、行列積などの基本的な行列計算のためのライブラリです。LAPACKなどの様々なプログラムから呼び出されます。行列計算は環境によって最適な書き方が違うので、行列計算の部分を差し替え可能なように分離しておきユーザが最適なBLASに差し替えて最高の性能を出す、ということを意図しています。

多くの環境で、デフォルトのBLASは単なるfor文で書かれたちっとも速くない「reference BLAS」と呼ばれるものになっています。自分の環境でも、reference BLASが入っていました。この場合例えばnumpyの行列計算は、全くCPU本来の性能を出せていません。

ubuntu 14.04では、atlasとopenblasという高速なBLASがパッケージで提供されています。今回は、それらをインストールしてnumpyの速度がどうなるか調べてみました。

インストールは次のように行いました。

sudo apt-get install libatlas-base-dev
sudo apt-get install libatlas-doc

sudo apt-get install libopenblas-base
sudo apt-get install libopenblas-dev

これで3種類のBLASがインストールされたことになりますが、これらは差し替えて使うことを意図しているのでshared libraryとして同一のファイル名になっています。Ubuntuでは、これらをsymbolic linkを張り直すことによって使い分ける仕組みがあります。

sudo update-alternatives --config libblas.so.3

と入力すると、

% sudo update-alternatives --config libblas.so.3
alternative libblas.so.3 (/usr/lib/libblas.so.3 を提供) には 3 個の選択肢があり>ます。

  選択肢    パス                                  優先度  状態
------------------------------------------------------------
* 0            /usr/lib/openblas-base/libblas.so.3      40        自動モード
  1            /usr/lib/atlas-base/atlas/libblas.so.3   35        手動モード
  2            /usr/lib/libblas/libblas.so.3            10        手動モード
  3            /usr/lib/openblas-base/libblas.so.3      40        手動モード

現在の選択 [*] を保持するには Enter、さもなければ選択肢の番号のキーを押してください:

のように表示され、どのBLASを使うか選択することができます。単にインストールした場合はUbuntuが定めた優先度に従って選ばれるようで、openblasが選択された状態でした。

これを切り替えてnumpyで5000x5000の逆行列を計算させてみて、速度がどうなるか計測してみました。環境は、core i7 3770T で、VMwareを使ってメモリを4G、CPUコアを4個与えた状態で行いました。

まずreference BLASで。

% python
>>> from numpy import *
>>> from time import *
>>> a=random.rand(5000,5000); t=time(); a=linalg.inv(a); time()-t
241.60200095176697

このように241秒かかりました。

ATLASに切り替えてみました。

% python
>>> from numpy import *
>>> from time import *
>>> a=random.rand(5000,5000); t=time(); a=linalg.inv(a); time()-t
32.59469819068909

32秒に高速化されました。しかし、topでCPUを見ると100%となっており、どうやら一つのコアしか使っていないようです。

最後にopenblasです。

% python
>>> from numpy import *
>>> from time import *
>>> a=random.rand(5000,5000); t=time(); a=linalg.inv(a); time()-t
6.014795780181885

6秒まで高速化されました。topのCPU表示は400%となり、4コア使って計算しているようです。

もともとテキサス大学の後藤先生という方が開発したgotoblas(ごとうブラス)という高速で有名なパッケージがあったのですが、後藤先生がIntelに移られて開発が中止してしまい、それを有志が引き継いだのがopenblasです。後藤先生は現在はIntelの商用BLASであるMath Kernel Libraryの開発をされているようです。

最後に、ATLASに関して少し注釈を述べておきます。
ATLASはAutomatically Tuned Linear Algebra Softwareの略で、ライブラリのmake時にそのCPUのコア数やキャッシュサイズなどを調査し、それに最適化したライブラリを生成する仕組みです。すなわち、本来バイナリパッケージで入れるには適さないもので、ATLASの真の性能を引き出すには計算を行うCPUでライブラリをmakeし直す必要があります。
また、ATLASはLAPACKの一部を高速化したものを含んでいて、ATLASをインストールするとATLAS版LAPACKがインストールされます。

% sudo update-alternatives --config liblapack.so.3
alternative liblapack.so.3 (/usr/lib/liblapack.so.3 を提供) には 2 個の選択肢があります。

  選択肢    パス                                    優先度  状態
------------------------------------------------------------
* 0            /usr/lib/lapack/liblapack.so.3             10        自動モード
  1            /usr/lib/atlas-base/atlas/liblapack.so.3   5         手動モード
  2            /usr/lib/lapack/liblapack.so.3             10        手動モード

現在の選択 [*] を保持するには Enter、さもなければ選択肢の番号のキーを押してください:

しかし、このように優先度が低くなっており、デフォルトでは使われていません。これがUbuntuの(というかDebianの)開発チームの見解ということでしょうか。

コメント（0件）

2015/06/06(土)kv-0.4.21

精度保証

kvライブラリを0.4.21にアップデートしました。

今回は、一応今までもひっそりと機能はあった、端点特異性を持つ関数の数値積分のコードを大幅に書き直し、機能も追加し、また詳細なwebドキュメントを書きました。端点特異性を持つ関数のうち、

除去可能特異点を持つもの。例えば $\int_0^1 \frac{\sin(x)}{x} dx$
ベキ型特異点を持つもので、x^pg(x)の形のもの。例えば $\int_0^1 \sqrt{x}\cos(x)dx$
ベキ型特異点を持つもので、f(x)^pの形のもの。例えば $\int_0^1 \sqrt{\sin(x)}dx$
ベキ型特異点を持つもので、f(x)^pg(x)の形のもの。例えば $\int_0^1 \sqrt{\sin(x)}\cos(x)dx$
log型特異点を持つもので、log(x)g(x)の形のもの。例えば $\int_0^1 \frac{\log(x)}{x+1}dx$
log型特異点を持つもので、log(f(x))の形のもの。例えば $\int_0^1 \log(\sin(x))dx$
log型特異点を持つもので、log(f(x))g(x)の形のもの。例えば $\int_0^1 \log(\sin(x))\cos(x)dx$

などを精度保証付きで計算できるようになりました。

その他、conj(共役複素数)を追加したり、全解探索における反復改良が稀にうまくいかなかったのを直したりしました。

最近少しずつユーザが増えてきて、最新のmacのclangで計算値が正しくないことがありそうとか、win7 32bit環境とVC++2012でddの表示がうまくいかないとか、報告を受けています。実機が用意できなかったりこちらでうまく再現しなかったりしてなかなか直せませんが、極力対応していきたいと考えていますので、こちらの環境(ubuntu 64bit)と異なる環境での動作報告、お待ちしております。

コメント（0件）