C言語からLAPACKのzgeevを呼んで複素行列を対角化
/* * SIZE*SIZE型の複素行列の固有値と固有ベクトルを計算 * (1 2i) * (i 1+i) */ #define SIZE 2 // 2*2型の行列 #include <stdio.h> #include <complex.h> int main(void) { char jobvl = 'N' ;// 左固有ベクトルは計算しない char jobvr = 'V' ;// 右固有ベクトルは計算する int n = SIZE ;// 対角化する正方行列のサイズ double _Complex A[SIZE*SIZE] ;// 対角化する行列 A[0]=1; A[2]=2*I; A[1]=I; A[3]=1+I; int lda = SIZE ;// 対角化する正方行列のサイズ double _Complex wr[SIZE] ;// 固有値が入る double _Complex vlDUMMY[1] ;// 左固有ベクトルが入る(不必要なので[1]) int ldvlDUMMY = 1 ;// vlDUMMY[1]の次元 double _Complex vr[SIZE*SIZE] ;// 右固有ベクトルが入る int ldvr = SIZE ;// 右固有ベクトルの本数 double _Complex work[2*SIZE] ;// 対角化する際に使用するメモリ int lwork = 2*SIZE ;// workの次元 double rwork[2*SIZE] ;// 2*SIZEで固定 int info ;// 成功すれば0 失敗すれば0以外を返す // LAPACKのzgeevサブルーチンを呼ぶ // 引数は全て参照渡し zgeev_(&jobvl, &jobvr, &n, A, &lda, wr, vlDUMMY, &ldvlDUMMY, vr, &ldvr, work, &lwork, rwork, &info); // 1番目の固有ベクトル : (vr[0] vr[1]) // 2番目の固有ベクトル : (vr[2] vr[3]) printf("1番目の固有値:%5.3lf %+5.3lf*I\n", creal(wr[0]), cimag(wr[0])); printf("1番目の固有ベクトル:(%5.3lf %+5.3lf*I %5.3lf %+5.3lf*I)\n", creal(vr[0]), cimag(vr[0]), creal(vr[1]), cimag(vr[1])); printf("2番目の固有値:%5.3lf %+5.3lf*I\n", creal(wr[1]), cimag(wr[1])); printf("2番目の固有ベクトル:(%5.3lf %+5.3lf*I %5.3lf %+5.3lf*I)\n", creal(vr[2]), cimag(vr[2]), creal(vr[3]), cimag(vr[3])); return 0; }
コンパイルと実行結果
$ gcc hoge.c -llapack -lblas -lm $ ./a.out 1番目の固有値:1.000 -1.000*I 1番目の固有ベクトル:(0.894 +0.000*I -0.447 -0.000*I) 2番目の固有値:1.000 +2.000*I 2番目の固有ベクトル:(0.707 +0.000*I 0.707 -0.000*I)
引数の詳細については、下手な解説を見るより zgeev.f の Arguments を読むのが一番わかりやすい。
CとFortranで行列の添字が異なる点への注意喚起
C言語
C言語の二次元配列 A[3][3]
とはつまり
A[0][0] A[0][1] A[0][2] A[1][0] A[1][1] A[1][2] A[2][0] ...
の順に並んだ一次元配列そのもの。 C言語の配列添字は0スタートだ。 従って、Cの標準的な二次元行列の表示は
( A[0][0] A[0][1] A[0][2] ) ( A[0] A[1] A[2] ) ( A[1][0] A[1][1] A[1][2] ) = ( A[3] A[4] A[5] ) ( A[2][0] A[2][1] A[2][2] ) ( A[6] A[7] A[8] )
Fortran
Fortranの二次元配列 A(3,3)
とはつまり
A(1,1) A(2,1) A(3,1) A(1,2) A(2,2) A(3,2) A(1,3) ...
の順に並んだ一次元配列そのもの。 Fortranの配列添字は1スタートだ。 従って、Fortranの標準的な二次元行列の表示は
( A(1,1) A(1,2) A(1,3) ) ( A(1) A(4) A(7) ) ( A(2,1) A(2,2) A(2,3) ) = ( A(2) A(5) A(8) ) ( A(3,1) A(3,2) A(3,3) ) ( A(3) A(6) A(9) )
数学
数学の配列添字は、まぁ大抵の場合1スタートだ。 二次元行列の数学的な表示は
a_{11} a_{12} a_{13} a_{21} a_{22} a_{23} a_{31} a_{32} a_{33}
結局どーやんの
僕達のやりたいことは、つまり
「数学的な二次元行列の{i,j}成分aij」を
「C言語のプログラム」上で
「Fortran式添字の変数A(i,j)に設置する」ことだ。
これは結局以下のようにすれば実現できると分かるだろう。
A[0+0*2] = a11; A[1+0*2] = a12; A[2+0*2] = a13; A[0+1*2] = a21; A[1+1*2] = a22; A[2+1*2] = a23; A[0+2*2] = a31; A[1+2*2] = a32; A[2+2*2] = a33;
訂正...コメントで間違いを指摘されました。ツメが甘かった。
/* C言語でFortran形式の行列を作る * * 数学的行列 * a11 a12 a13 * a21 a22 a23 * a31 a32 a33 * * を一次元配列 A[0] ... A[8] に入れる。 * A[(i-1)+(j-1)*3] に aij を入れれば良い。 */ double _Complex A[3*3]; A[0+0*3] = a11; A[0+1*3] = a12; A[0+2*3] = a13; A[1+0*3] = a21; A[1+1*3] = a22; A[1+2*3] = a23; A[2+0*3] = a31; A[2+1*3] = a32; A[2+2*3] = a33;
2017/01/20 追記
zheevの使い方を書きました。 auewe.hatenablog.com