2015/09/10(木)応用数理学会年会
2015/08/06(木)kv-0.4.23
今回は、allsol.hppにある非線形方程式の全解を探索するプログラムに、[0,∞]や[-∞,∞]などの無限大を含む領域を探索区間として指定できるようにしました。もっとも、必ず全ての解を求めて停止するわけでは無く、「求まる問題なら求まる」といったレベルです。ある問題を解いていてたまたま必要になったので実装したのですが、まだまだ改良が必要と思います。こうして文章にすると簡単そうですが、実は結構苦戦しました。無限大は難しい。
その他、vleq.hppの仕様を少し変更したり、区間のmidを精密にしたりしました。ベキ級数演算のpdfのドキュメントにも少し加筆しました。
2015/06/29(月)kv-0.4.22
今回は、defint_power_r, defint_log_r, defint_power_autostep_r, defint_log_autostep_rの4つの端点特異性を持つ積分の関数のbug fixです。これら右側の端点特異性を扱う関数は、時間を反転させて左側versionに任せているのですが、その時間反転部分にバグがありました。これらの関数を使う例題を作っていなかったので、気付きませんでした。
その他、細かいところをいくつか修正しています。
2015/06/28(日)octave+openblas+intlabで丸め変更を検証
行列積はn3回の乗算とn2(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上で実行してみました。なお、数値は行列サイズで、このプログラムはn4に比例した時間がかかるのであまり大きな行列ではテストできません。実行結果は、次のようになりました。| 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 |
# 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
ところで、INTLAB公式のoctaveページからダウンロードしたoctave-3.8.2-installer.exe、Aviraさんが
みたいに反応するんだが、大丈夫だろうか (自分は仮想環境で実行してすぐに消しました)。Aviraさんはときどき過剰に反応するからそのせいだと信じよう。
2015/06/24(水)INTLAB version 9をubuntu上のoctaveで使う
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
実際に区間演算をしてみます。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の使い方はよく分かっていませんが、お楽しみ下さい。