// dsbev.cc
// 実帯行列の固有値方程式を解く

#include <stdio.h>

extern "C" {
  void dsbev_ ( const char& JOBZ, const char& UPLO,
		const int& N, const int& KD,
		double** AB, const int& LDAB,
		double* W,
		double** Z, const int& LDZ,
		double* WORK,
		int& INFO, int JOBZlen, int UPLOlen );
};

// 実対称3重帯行列 H の固有値方程式 H x = E x を解く簡易関数
// Input: H[2][N],  Output: U[N][N], E[N]
// H[0][N] に対角成分、H[1][N] 下副対角成分
// U^\dagger H U が対角行列となる
//
template <int N> int dsbev3( double H[2][N], double U[N][N], double E[N] )
{
  int i, j, info;
  static double work[4*N];
  static double AB[N][N];

  for( j=0; j<N; j++ ){
    for( i=0; i<2; i++ ){
      AB[j][i] = H[i][j];
    }
  }

  dsbev_( 'V', 'L', N, 1, (double**)AB, N, E, (double**)U, N, work, info, 1, 1 );

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

  return info;
}

// 実対称5重帯行列 H の固有値方程式 H x = E x を解く簡易関数
// Input: H[3][N],  Output: U[N][N], E[N]
// H[0][N] に対角成分、H[1,2][N] 下副対角成分
// U^\dagger H U が対角行列となる
//
template <int N> int dsbev5( double H[3][N], double U[N][N], double E[N] )
{
  int i, j, info;
  static double work[4*N];
  static double AB[N][N];

  for( j=0; j<N; j++ ){
    for( i=0; i<3; i++ ){
      AB[j][i] = H[i][j];
    }
  }

  dsbev_( 'V', 'L', N, 2, (double**)AB, N, E, (double**)U, N, work, info, 1, 1 );

  for( i=0; i<N; i++ ){
    for( j=i+1; j<N; j++ ){
      double 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;
  double H[3][N], U[N][N], E[N];

  for( j=0; j<N; j++ ){
    H[0][j] = 0.0;
    H[1][j] = -1.0;
    H[2][j] = -0.5;
  }

  info = dsbev5( H, 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 ", U[i][j] );
    }
    printf("\n");
  }

  return 0;
}