CLAPACK使ってみた

取り敢えずサンプルを動かすのには成功しました。参考にしたのは、

http://www.k.mei.titech.ac.jp/~stamura/NumericalComputation-Tips.html#compile-LAPACK
http://www.cazoo.jp/calcsoft/library/lapack/clapack/usage/dgeev/index.html

で、上の方はコンパイルのときのオプションのつけ方。下の方はサンプルのコードがある感じ。俺はコンパイルは、

$gcc hoge.c -o hoge -lm -llapack -lblas -lg2c

にした。サンプルコードは丸ごとコピペで、

#include <stdio.h>
#define N 3
double A[N*N];
double wr[N];
double wi[N];
double vl[1*N];
double vr[N*N];
double work[4*N];
int main(void)
{
  static long int i;
  char jobvl = 'N';
  char jobvr ='V';
  static long int n=N,lda=N,ldvl=1,ldvr=N,lwork=4*N,info;

  A[0]=6.;A[1]=-1.;A[2]=5.;
  A[3]=-3.;A[4]=2.;A[5]=-3.;
  A[6]=-7.;A[7]=1.;A[8]=-6.;

  dgeev_( &jobvl, &jobvr, &n, A, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info);

  for(i=0;i<N;++i) printf("%lf\n",wr[i]);

  return 0;
}

こんな感じ。これが行列A[i][j]の固有値を求めるコードらしい。何か右固有値とか左固有値とかあるらしいです。右固有値が、普段固有値と言ってる、
  A\vec{x}=\alpha\vec{x}
のαで、左固有値っていうのは、上の式の転置をとった、
  \vec{x}^tA^t=\beta\vec{x}^t
のβのことらしいです。量子力学とか連続体力学やる人は転置とるのを考えたりするのかね?連続体力学はこういうのが沢山出てきて目がチカチカした記憶がある。全然分からんかったよ。俺には流体が限界ですよ。そして量子力学ではHermiteがどうとかやったなあ。懐かしいなあ。っつーかHermite行列はかなりあり難い感じです。対称行列とかね。これらの性質があるからまだ量子力学流体力学もなんとかできてる感じがする。っつーかまあそもそも世の中が非対称だったらそれはそれで困るわけだけど*1
そして俺がやるのはこれなので、これから使い方を勉強して、今組んでるコードに実装します。
そいでさー、なんつーかさあ、世の中にはLAPACKについて言及してあるサイトは糞みたいに沢山あるんだけど、具体的にどんなヘッダをインクルードして、どうやってサブルーチンを呼び出して、コンパイルのときにどういうオプションつければ良いみたいな情報が殆どないんですよね。そんな固有値問題の説明解かされてもマジで困るんですよね。
恐らく普段使いこなしてる人にとってはパッケージの使い方なんざ自明なことなので書かなくても分かるんだろって感じなんだろうけど、こっちゃ初心者だし、Webに公開するなら初心者が見ることを念頭に置かないとそれが価値のある情報であると判断されるかってすげー微妙だと思うんですよね。まあ俺も人のことはいえませんが。
っつーことで色んなものへの不信感が募って妙に疲れた。っつーか仕様だけ書かれてもわかんねーっつーの。こっちは馬鹿なんだっつーの。導入から書けやヴォケがって感じですね。往々にして世の中の理系の自称とっても良くできる人は他人が何が分からないかとかいうことへの配慮が欠けてる気がします。まあ俺も人のこと言えないんですけど。
因みにCLAPACKについては、俺が参考にした

http://www.cazoo.jp/calcsoft/library/lapack/clapack/index.html

がとっても良くまとまってて分かり易いです。ここまでやるのは相当なことかなあと。こういうのがあるとマジ助かる。

*1:これは剪断応力とか、変形速度テンソルが対称っていう議論についてです。量子力学で、状態量が非Hermiteだったらどうなるんだっけか?すっかり忘れた。っつーかCANNAは「剪断」て出力できるのが良いですね。