// zhegv.cc
// 複素エルミート行列の一般化固有値方程式を解く

#include <stdio.h>
#include <complex>
typedef std::complex<double> Complex;

extern "C" {
  void zhegv_ ( const int& ITYPE, const char& JOBZ, const char& UPLO,
		const int& N,
		Complex** A, const int& LDA,
		Complex** B, const int& LDB,
		double* W, Complex* WORK, const int& LWORK,
		double* RWORK,
		int& INFO, int JOBZlen, int UPLOlen );
};

// 複素エルミート行列 H,S の一般化固有値方程式 H x = E S x を解く簡易関数
// Input: H[N][N], S[N][N],  Output: U[N][N], E[N]
// ただし、H[i][j<=i] S[i][j<=i] の下三角しか参照されない
// U^\dagger H U が対角行列、U^\dagger S U が単位行列
//
template <int N> int zhegv( Complex H[N][N], Complex S[N][N], Complex U[N][N], double E[N] )
{
  int i, j, info;
  static Complex cwork[4*N];
  static double  rwork[4*N];
  static Complex B[N][N];

  for( i=0; i<N; i++ ){
    for( j=0; j<N; j++ ){
      U[i][j] = H[j][i];
      B[i][j] = S[j][i];
    }
  }

  zhegv_( 1, 'V', 'L', N, (Complex**)U, N, (Complex**)B, N, E, cwork, 4*N, rwork, info, 1, 1 );

  for( i=0; i<N; i++ ){
    for( j=i+1; j<N; j++ ){
      Complex temp;
      temp = U[i][j];
      U[i][j] = U[j][i];
      U[j][i] = temp;
    }
  }

  return info;
}

int main(void)
{
  const int N = 4;
  int i, j, info;

  Complex H[N][N], S[N][N], U[N][N];
  double E[N];

  for( i=0; i<N; i++ ){
    for( j=0; j<N; j++ ){
      if( i == j+1 || i == j-1 ){
	H[i][j] = -1.0;
      }else{
	H[i][j] = 0.0;
      }
      if( i == j ){
	S[i][j] = 1.0;
      }else if( i == j+1 || i == j-1 ){
	S[i][j] = 0.5;
      }else{
	S[i][j] = 0.0;
      }
    }
  }

  info = zhegv( H, S, U, E );

  printf("# info=%d.\n", info );

  printf("# Eigen values.\n");
  for( j=0; j<N; j++ ){
    printf("%+f ", E[j] );
  }
  printf("\n");

  printf("# Eigen vectors.\n");
  for( i=0; i<N; i++ ){
    for( j=0; j<N; j++ ){
      printf("%+f%+f ", real(U[i][j]), imag(U[i][j]) );
    }
    printf("\n");
  }

  return 0;
}