2015年の数値計算日記
こけつまろびつ

桂田祐史

Date: 2015年9月1日, 2018年3月1日

去年, 先月, 今月, 来月

12月28日(月曜)

私事が忙しいとはいえ、日々の授業の準備を除くとあまり数学をしていない。一応冬休みになって、それでも私事の合間を縫ってやったのは、自由な縁を持つ円形板のクラドニ図形。 19世紀に出来ていたことだけど、 Courant-Hilbert では大まかなことしか書いてなくて、グチャグチャした (でもまあまあ詳しい) Rayleigh を読み出した。式の計算は良いけれど (チェック出来た)、計算出発の前提である、固有関数が

$\displaystyle \cos(m\theta+\phi)\left[J_m(kr)+\lambda I_m(kr)\right]$

という形をしているというのが良く分からなかった。常微分方程式でないので、 $\Laplacian^2-k^4= \left(\Laplacian^2-k^2\right) \left(\Laplacian^2+k^2\right)$ という因数分解で問題を分解して良いものか。いろいろ考えて、金子晃先生の本で学んだ考え方で乗り切れるのかも、と思いついた。詰めるのにはもう一努力が必要か。
ともあれ、Rayleigh が何を主張しているのかは明確に分かったので、それを計算するコードを書けるようになった。紆余曲折あったけれど、出来てしまえば大筋は棒の固有値問題の計算と同じ。
クラド二の法則は、円形板で、クラド二が $m\vert n$ という記号を割り当てた振動モードの周波数についての近似法則であり、

$\displaystyle f_{mn}\simeq C(m+2n)^p, \quad p\simeq 2$

というもの。ところで $f_{mn}=C k_{mn}^2$ なのだけど、 ``ありがちな''

$\displaystyle k_{m+1,n}-k_{m,n}\to 1,\quad k_{m,n+1}-k_{mn}\to \pi$

が見えてきて (多分証明出来る)、「あれれ」となる。のは $\pi$ なんじゃないかな。これはネタになるのかも。
　そうそう。 RIMS 研究集会「偏微分方程式の逆問題とその応用の新展開」で 1/29(金曜) 13:45 に藤原さんが ``Multiple-Precision Arithmetic on MATLAB and its Application to the Real Inverse Laplace Transform'' という講演をする。聞きに行く？来年の年会は筑波なのでそれにお金は必要ない。
　そうそう。ソフトウェアセミナー：FreeFem++による有限要素プログラミング -- 中級編 --

11月20日(金曜)

昨日は総合数理概論の準備を割と早くスタートして、色々手を入れた。固有振動のプログラムを、 C+GLSC から Mathematica に切り替えた。
c=1; L=1; u[n_, x_, t_] := Sin[n Pi x/L] Cos[c n Pi t/L] Animate[Table[ Plot[{u[n, x, 0], u[n, x, t]}, {x, 0, L}, PlotRange -> {-1, 1}], {n, 1, 3}], {t, 0, 2L/c, 0.01}]

11月11日(木曜)

久しぶりに「単振り子の話」に加筆した。 http://nalab.mind.meiji.ac.jp/~mk/labo/text/furiko/node12.html

10月21日(水曜)

夜中に考えた。近似解析関数が得られることと、その極や零点を得ることとの間には差があると。それは言わないといけない。
そのうちに考えた。それは本当だろうか？ Rouché の定理を考えた。ああいうのも周回積分だ。
関数論は周回積分で構築されている。解析的曲線に沿う周回積分は台形則で効率的に近似出来る。
そういうのを使えば色々出来るのではないか？
例えば、微分も Cauchy の積分公式によれば周回積分表示できる。
もしも精度保証付き数値計算が出来れば凄いかも。
数値微分とどっちが良いか考えた。あれ、円周上に二点を取ると普通の中点則なのか？その観点から調べてみると面白いかも。
CDSM も合流できるのでは？
もしかして、 $\dsp\int_{\left\vert z-c\right\vert=R}z^n\;\Dz=0$ ( $n\ne -1$ ) は、普通に台形則近似して 0 にならんか？
そうそう。楕円関数の数値計算を考えておこう。
風上差分の質問を受けた。そういうの勉強したいな。
移流方程式の初期値境界値問題。

    $\displaystyle u_t=c u_x$ ( $x\in(0,\ell)$ , ) $\displaystyle ,$

    $\displaystyle u(0,t)=a$ ( ) $\displaystyle ,$

    $\displaystyle u(x,0)=u_0(x)$ ( $x\in[0,\ell)$ ) $\displaystyle .$

一瞬 well-posed かどうか気になったけれど。なるほど、これは大丈夫そうだ。

$\displaystyle u(x,t)= \left\{ \begin{array}[t]{ll} u_0(x-ct)& \text{($x-ct>0$)} \\ a & \text{($x-ct<0$)} \end{array} \right.$

だろうか？
I(M1)君にお説教。それで結合振動子の話になる。本を少し注文。メトロノームの同期の話を調べたら、龍谷の池田先生の学生の論文とかヒット。面白い。

10月20日(火曜)

1ヶ月飛んだ。なぜだろう。何だかんだ刺激がないということだ。
N先生から共同研究のお誘いがあって、それで打ち合わせ。飛んでる話だ。
複素変数の常微分方程式はよくある。
Cauchy-Kowalevski の定理があるな。
双曲型方程式の初期値問題の一意性は何か解析性と関係があるか？特異性の伝播とか。
楕円型偏微分方程式の解は一般に実解析的なのか？
楕円型方程式の一意接続定理というのがあった。どんなの？

9月19日(土曜)

なぜか MacPorts で Octave がうまく make 出来ない。というのと、配布されているバイナリーでも、 MATLAB で動いているプログラムがきちんと動かない。ずっしりと重い処理をさせているせいなのかもしれないけれど。
Octave とは気長に付き合うことにするけれど、 MATLAB+INTLAB についても調べよう。
MATLAB+INTLAB については、怪しい噂も聞く(読む)けれど、具体的な事実を書いていない噂だけで判断しても仕方がないので、まずはやってみましょう。
MATLAB R2015a をインストールした。一度 MATLAB を起動すると、 ドキュメント/MATLAB というディレクトリィが出来る。

INTLAB をインストールしよう。まずターミナルから

unzip INTLAB.zip
cp -pr Intlab_V9 ~/Documents/MATLAB

MATLAB のコマンド・ウィンドウで

cd Intlab_V9
startintlab

次のようなメッセージが表示される。何気に大事そうなことが書いてあるので、良く読むこと。

VVVVVVVVVVVVVVVV Installation of INTLAB	VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
 
This is the first call of INTLAB, and the best way to switch the 
rounding mode will be detected. It is a very rigorous check and takes
a short moment. This is done only once with the very first call of
		INTLAB.
 
To find the best routine, various routines are checked. 
Unfortunately, in very rare cases it might happen that a hard error
occurs. This should never happen, but it happes in some Matlab versions
and/or operating systems, and it is beyond my control. 
 
In such a rare case of a hard error, please restart INTLAB.
The circumstances will be stored, and this switching routine will not
be called again. Sorry for possible inconveniences. 
 
Press Enter to continue.

使うたびに PATH を気にしたくないので、設定ファイルに PATH を加えるように書いておく。

~/Documents/MATLAB/startup.m に書き加える

addpath('~/Documents/MATLAB/Intlab_V9')

こうしておけば、MATLAB から

>> startintlab

で INTLAB が起動する。

丸めの制御が出来るかどうかは、お金を支払う前に分かるようにすべきでは？という気がしてきた。

18日(金曜)

本当に新しい MATLAB では INTLAB がきちんと動かないのだろうか？もしそうならば、そのことを明記しないで有料 (50ユーロ) で配布しているのは、ちょっと問題があるような気がする。
INTLAB 自体に自己チェックの機能があるはずで、新しい MATLAB でもインストールして試してみるか。
精度保証付き数値計算の資料ページに少し手を入れる。
文献みたいのも入手すべきか。まずは中尾・山本の電子化をするとか (あ、アマゾンで 5000 円になっている, まあ持っているから良いか)。中尾・渡部も買う (電子化用、本棚用に2冊)。
大石先生は最近は本は出していないのかな。発展中だしね。まあ、線形計算の方は中尾理論ほど敷居が高くないので、「論文を読め」で良いのかな。

17日(木曜)

ある時期、線形計算の精度保証付き数値計算について、研究の現状の理解をしないと、と考えたことがあったけれど、忙しさに紛れて有耶無耶になってしまった。テキストのない状態で、今から追いかけるのは大変のような気がするが、何をするにも線形計算は無視できない。逃げることは出来ない。
というわけで論文の収集をすることに。最近の発表のレファレンスを見て入手すべき資料を拾い上げるのかな。何か勉強ノートみたいのを作るのか。最初は非公開で。学生が誰か卒研でやってくれると深まるかも。

15日(火曜)

MacPorts について色々調べた。 octave に openblas という variant を付け加えられれば良い。 openblas にも consistentfpcsr という variant を加える?
岡本久さんの総合講演のアブストラクトを読んだ。数学にとっても数値計算の意味ということを深く考えさせられた。電子化して学生に読んでもらおう。

14日(月曜)

OpenBLAS 昨日のまとめ。とりあえず自前で作る。

tar xzf OpenBLAS-0.2.14.tar.gz
cd OpenBLAS-0.2.14
vi Makefile.rule

#CONSISTENT_FPCSR = 1 をアンコメントする。

make>&make.log
sudo make PREFIX=/usr/local install

Octave 昨日のまとめ。「大変だ。」じゃなくて、
- MacPorts では、 BLAS について、 atlas 使うか、バギーな vecLib 使うかの二択になっている。 OpenBLAS を使いたい。途中まで vecLib にしておいて、最後に -L某 -lopenblas とすると良いのでは？
- sudo port install octave +accelerate +gcc5 +glgui で作業を始める。
- port dir octave とすると、 /opt/local/var/macports/sources/rsync.macports.org/release/tarballs/ports/math/octave と答が来る。
  cd /opt/local/var/macports/sources/rsync.macports.org/release/tarballs/ports/ cd math/octave/work/octave-3.8.2
  でコンパイルする作業ディレクトリに移れる。
作業後に消えないようにしたいが、そうでなければ、Building と出たら C-z で中断するのだろうか？必要ならばコピーしてから fg とする。
sudo port extract octave +accelerate +gcc5 +glgui という手もあるらしい。
まず作ってから BLAS の設定を変えて make してみる。自分で +openblas というのを作れると良いね。
相変わらず TeX 関係で引っかかるので、そのタイミングで make clean するのかな。
もしかして -with-blas=-L/usr/local/lib -lopenblas とすれば良いのかな。

13日(日曜)

現在、Mac を使っていてなかなか快適なのだけど、精度保証付き数値計算や、多倍長計算をするための環境は未整備である。そこら辺をコツコツ作業していこう。
どういうようにしたいか、いくつか書いてみる。
- Octave を使えるようにしよう。 -- 十数年前は MATLAB+INTLAB がそびえ立っていたが、最近 MATLAB 側の都合でうまく動かないらしい。また動くようになることを期待するが、 Octave で INTLAB が使えるそうなので、まずは Octave をしてみよう。
- 多倍長で精度保証付き計算をする必要がありそうだが、 Octave や MATLAB で出来るようになるのは時間がかかりそうで、 C++ や Fortran で使えるようなライブラリィの導入が必要だ。
- 以前、多倍長計算用のクラス・ライブラリィはいくつか試してみたが、区間演算が出来るものはあるだろうか？
- kv は作者の哲学も良さそうで、試してみよう。
- とすると Boost がいるのか？
- 多分 BLAS は必須だろう。ATLAS は何となく良い印象がないので、 OpenBLAS なんですかね。
- それ自体を使うかどうか良く分からないが、LAPACK を入れる。 BLAS のテストもし易いだろうし。そのために Fortran を入れる必要が出て来るが、大して邪魔にもならないだろう。
- もしかすると Mac でなくて Linux なのかも。
まずは、LAPACK, BLAS かな。どちらかというと、 LAPACK よりも BLAS に重点がある。精度保証付き数値計算に使える、高い効率の BLAS を、というのが目標である。
Cコンパイラー、Fortran コンパイラーについて。 Xcode はかつては GCC を採用していたが、最近ではそれをやめて、LLVM が採用されている。 GCC と LLVM, どちらを選ぶべきか実は良く分かっていないが、とりあえず GCC (GNU Compiler Collection) を採用してみる。現在、元気に開発が進んでいる kv などで、メインに採用されていそうだから、という安直な理由である。
gcc, gfortran を自分で一から make する気になれない。 HPC Mac OS X から、 gfortran, gcc のバイナリーを入れるという手もあるが、年に1回しかリリースされないようだし、今だと Mavericks にインストールできるのは、 gcc-4.9-bin.tar.gz, gfortran-4.9-bin.tar.gz になる。gcc-5.1 は Yosemite のみだ。
これに対して、MacPorts ならば、gcc-5.2 がインストール出来る。まあ、MacPorts で行こう。
sudo port install gcc5
コマンド名は gcc-mp-5, gfortran-mp-5 になるようだ。すでに自分のマシンにはインストールされていたようで、 gcc5 だけで足りるのかは良く分からない。
BLAS が必要だ。OS X に入っているのは評判が悪いので、OpenBLAS か。 OpenBLAS によると、 User Manual と Installation Guide を読め、とある。
If your CPU is Sandy Bridge, please use Clang version 3.1 and above. The Clang 3.0 will generate the wrong AVX binary code of OpenBLAS.
えー？Core i7 というのはそれに入るのか？ Clang のバージョンばかり気にしているようだけど、 GCC は使うなということなのかな。
http://github.com/xianyi/OpenBLAS/archive/v0.2.14.tar.gzを入手して、
tar xzf OpenBLAS-0.2.14.tar.gz cd OpenBLAS-0.2.14 make>&make.log
エラーになった。あれれ。ログを見る。 g95 を使っているぞ。何かおかしい。 g95 が /opt/local/bin にあったけれど、タイムスタンプ 2010 年だって。怪しい。フル過ぎる (/opt/local/bin にあるもので一番古い)。
g95 を削除してから、make したら、またエラー。何となく make clean してから、再度 make したら通った。つまり古い g95 がなければすんなり通ったはずだ。 gfortran でなく、g95 が優先されたのは謎だ。
どこにインストールするのか迷ったけれど、とりあえず /usr/local に。
sudo make PREFIX=/usr/local install
反省してみる。 octave+openblas+intlabで丸め変更を検証にあるように OpenBLAS の Makefile.rule
# CONSISTENT_FPCSR = 1
をコメント・アウトした。こうしても make できた。
以上は clang でコンパイルされていた。 OpenBLASは gcc で make すべきでは？この線を検討してみる。
make CC=gcc-mp-5 FC=gfortran-mp-5
とすると
Cannot find plugin 'liblto_plugin.so'
というエラーで止まる。
ln -s /opt/local/bin/ar gcc-ar ln -s /usr/bin/ar gcc-ar
試しに MacPorts で openblas をインストールしてみたら、 LLVM 3.4 がインストールされた。うーんんん。システムの clang が LLVM 3.5 を使っているというのに。
DragonEgg - Using LLVM as a GCC backend ということか。
色々なものの集合体ということは、一つの目的を持っている人からすると、段々ずれていく可能性があるということかな。
Makefile.rule をきちんと書き換えろ、という説も。そういう人が書いた Makefile.rule を見ても、デフォールトで良いはずのところしか書いていない。
まあ、clang でコンパイルして、ちゃんと動けば良いわけだよな。 OpenBLAS は Microsoft のコンパイラーでもきちんと動くらしいし、 GCC オンリーというわけではない。 CONSISTENT_FPCSR = 1 が有効になれば文句はないのかも。
port dir openblas として、 /opt/local/var/macports/sources/rsync.macports.org/release/tarballs/ports/math/OpenBLAS/work/OpenBLAS-0.2.14 を知る。そこの Makefile.rule に CONSISTENT_FPCSR = 1 を書き足した。

まとめておく。

素直に make (ただし CONSISTENT_FPCSR = 1 だけはしておく。)

tar xzf OpenBLAS-0.2.14.tar.gz
cd OpenBLAS-0.2.14
vi Makefile.rule
make>&make.log
sudo make PREFIX=/usr/local install

/usr/local にインストールされる。

MacPorts でインストール

sudo port install openblas
 C-z (中断)
pushd /opt/local/var/macports/sources/rsync.macports.org/release/tarballs/
cd ports/math/OpenBLAS/work/OpenBLAS-0.2.14
vi Makefile.rule

(CONSISTENT_FPCSR = 1 を有効にする。)

popd
fg

/opt/local にインストールされる。これも clang でコンパイルされている。

gcc, gfortran が検索できるようにしておいて、 CC=, FC= で指定する。 (そうしないと、clang を選択するね。Makefile.rule に書くべきか？)
tar xzf OpenBLAS-0.2.14.tar.gz cd OpenBLAS-0.2.14 vi Makefile.rule make CC=gcc FC=gfortran>&make.log sudo make PREFIX=/usr/local install
やはり sgemv_t_microk_sandy-4.c が扱えないと言う。

気分転換に LAPACK を。 /opt/local/lib に liblapack.a があった。日付から見て arpack で必要になったからインストールされたらしい。まあ /usr/local/lib にインストールするかな。

LAPACKから lapack-3.5.0.tgz を入手する。

tar xzf lapack-3.5.0.tgz cd lapack-3.5.0 cp INSTALL/make.inc.gfortran make.inc vi make.inc
(特に書き換えるべきところはない？以前は PLATFORM を LINUX から DARWIN にしたりしたけど。)
make blaslib variants>&make-blaslib-variants.log make >&make.log

コンパイルに警告１つ出ない。テストも全部通った。

sudo cp -p liblapack.a /usr/local/lib/ sudo cp -p libtmglib.a /usr/local/lib/

LAPACK の C インターフェイスである LAPACKE は、 LAPACK に取り込まれたそうだ。

cd lapack-3.5.0
make lapackelib

それから

sudo cp -p liblapacke.a librefblas.a /usr/local/lib
sudo cp -p lapacke/include/*.h /usr/local/include

さて、テストだ。 test-lapacke.c

gfortran -I/usr/local/include test-lapacke.c \
 -L/usr/local/lib -lmatrix -llapacke -llapack -lopenblas

% ./a.out
A=
  11   12   13    0    0    0    0    0
  21   22   23   24    0    0    0    0
  31   32   33   34   35    0    0    0
   0   42   43   44   45   46    0    0
   0    0   53   54   55   56   57    0
   0    0    0   64   65   66   67   68
   0    0    0    0   75   76   77   78
   0    0    0    0    0   86   87   88
b=
74 230 505 890 1385 1990 1994 1829

call dgbsv
 successful
 return from dgbsv

計算して得た解=
1 2 3 4 5 6 7 8
 結果が 1,2,3,.. となっていたら成功
%

今のところ、ボロなく動いている模様。めでたい。
MacPorts に BLAS/LAPACK がないのは、OS X のを信じているの？？
さて、Octave やってみるか。ATLAS は使わないで。とりあえず +accelerate でやってみる？何もしないでインストールを始めて、
port variants octave
sudo port install octave +accelerate +gcc5
ちなみに /opt/local/var/macports/sources/rsync.macports.org/release/tarballs/ports/math/octave/work/octave-3.8.2 で作業をしている。ここをちょっとコピーして後でいじることに。
Octave は BLAS をダイナミック・リンクするので、コンパイルし直す必要がないという説もある。

あ、ひっかかった。

cd /opt/local/var/macports/build/
cd _opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_octave
cd octave/work/octave-3.8.2
sudo make
sudo port install octave +accelerate +gcc5

実行してみる。なんか変だな。相変わらずダメだ。

sudo port install octave +accelerate +gcc5 +x11

バイナリーで配布されている Octave
otools -L /usr/local/octave/3.8.0/bin/octave-3.8.0
BLAS なんて出て来ないな。静的にリンクされているのだろうか。

sudo make check

  libinterp/corefcn/dlmread.cc-tst ....................... PASS   12/20	FAIL 8

  libinterp/corefcn/sparse-xpow.cc-tst ................... PASS    1/2	FAIL 1

  libinterp/corefcn/str2double.cc-tst .................... PASS   28/31	FAIL 3

  scripts/io/importdata.m ................................ PASS   22/23	FAIL 1

  scripts/sparse/svds.m .................................. PASS    5/6	FAIL 1

  scripts/specfun/realpow.m .............................. PASS    4/5  FAIL 1

  io.tst ................................................. PASS   88/91	FAIL 3

The HDF5 header files used to compile this application do not match
the version used by the HDF5 library to which this application is linked.
Data corruption or segmentation faults may occur if the application continues.
This can happen when an application was compiled by one version of HDF5 but
linked with a different version of static or shared HDF5 library.
You should recompile the application or check your shared library related
settings such as 'LD_LIBRARY_PATH'.
You can, at your own risk, disable this warning by setting the environment
variable 'HDF5_DISABLE_VERSION_CHECK' to a value of '1'.
Setting it to 2 or higher will suppress the warning messages totally.
Headers are 1.8.15, library is 1.8.14

General Information:
-------------------
		   HDF5 Version: 1.8.14
		  Configured on: Tue Jan 13 07:47:08 CET 2015
		  Configured by: hecht@MBP-FH.local
		 Configure mode: production
		    Host system: i386-apple-darwin14.0.0
	      Uname information: Darwin MBP-FH.local 14.0.0 Darwin
		Kernel Version 14.0.0: Fri Sep 19 00:26:44 PDT 2014;
		root:xnu-2782.1.97~2/RELEASE_X86_64 x86_64
		       Byte sex: little-endian
		      Libraries: static, shared
	     Installation point: /usr/local

2015年2月19日って何？もっともリンクは6月21日だ。GLSC? あ、FreeFem++ かな。
UMFPACK, ARPACK, rl_set_keyboard_input_timeout CXSparse CHOLMOD COLAMD CAMD AMD FFTW3F
-with-hdf5-includedir=/opt/local/include -with-fftw3-includedir= -with-fftw3-libdir
./configure -prefix=/opt/local -includedir=/opt/local/include -libdir=/opt/local/lib

sudo port install octave +accelerate +gcc5 +glgui

12日(日曜)

思い立って、今日から「数値計算日記」を書くことにした。
もともとコンピューター上で日記を書いていたのだが、プライベートなこと (自分や大学に関すること) は除いて、数値計算に関する勉強や環境整備などに関することだけ別立てにしよう、ということである。
そうするつもりになった理由は二つある。
1. いくつかの事項 (1つには精度保証付き数値計算) について、まとまった勉強をしよう。
2. ゼミの学生に数値計算に関する情報を何らかの形で伝えたい。
精度保証付き数値計算はその道の専門家でないと分かりにくくなってしまったな、これからは野次馬するだけに止めておくものか、と感じたこともあったけれど、どうもとんでもない勘違いをしていたようだ。すべての種類の数値計算について、可能ならば精度保証をする方向に進むような気がしてきた。自分が生きている間、ボケない間だけでもかなり進展するのではないだろうか。
というわけで気持ちを改めて勉強を始めることにした。

この文書について...

桂田祐史
2018-03-01

		$\displaystyle u_t=c u_x$ ( $x\in(0,\ell)$ , ) $\displaystyle ,$
		$\displaystyle u(0,t)=a$ ( ) $\displaystyle ,$
		$\displaystyle u(x,0)=u_0(x)$ ( $x\in[0,\ell)$ ) $\displaystyle .$

2015年の数値計算日記 こけつまろびつ

2015年の数値計算日記
こけつまろびつ