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