やったー、やったよままん、自分でコード組んでCLAPACK動かせたよ。固有値求めるところまではうまくいったんだけど、その後の固有ベクトルを求めるところが上手く行かなくて、ずーっと試行錯誤してて、最終的に作業領域の大きさの指定の不整合で上手くいってないんじゃないかと思うようになり、結局それを直してみたら上手くいったよ。
っつーか俺はネットから落としてきたサンプルコードをそのままコンパイルして実行したんだけど、上手く行かなかったんだよね。で、結局そのサンプルは微妙に間違えててみたいな感じで。っつーか配列の切り方とかマジで重要です。あんなのちょっと間違えただけででっかいミスになるなんて阿呆臭くてやってらんないっすよ。
っつーことで、俺の組んだコードを公開してみます。
の固有値と固有ベクトルを求めるプログラムです。
#include <stdio.h> #include <math.h> #include <complex.h> #define N 2 int main () { long int dim = N, lda = N, ldvl = 1, ldvr = N, lwork = 2*N, info; double _Complex A[N*N], EValue[N], EVector[N*N]; double _Complex work[2*N], rwork[2*N]; int i; A[0] = 1; A[1] = 1+I; A[2] = 1-I; A[3] = 1; zgeev_("N", "V", &dim, A, &lda, EValue, NULL, &ldvl, EVector, &ldvr, work, &lwork, rwork, &info); printf("info=%ld\n", info); for (i = 0; i < N; i++) { printf ("EigenValue%1d=%4lf+%4lfI\n", i+1, creal(EValue[i]), cimag(EValue[i])); } printf ("Eigen Vectors are shown bellow\n"); for (i = 0; i < N*N; i++) { printf ("%lf+%lfI\n", creal(EVector[i]), cimag(EVector[i])); } return 0; }
もとは、
http://www.cmpt.phys.tohoku.ac.jp/~ihasato/HowToUseCLapack/zgeev0.cpp
これを参考にしてるんだけど、これが微妙に間違えてたので、自分なりに直してみました。そしたら動いたよ。マジで感動した。いやープログラミング、嫌いだけど、こういうのがあるから本当に嫌いにはなれないんだよなあ。
あとはまあCLAPACKについては資料が少ないんで、困ってる人もいるかも知れないっつーことで、ここに載せることにしました。
そして今日の日記は糞長いですね。