// zgbsv.cc
// 複素帯行列の線形方程式を解く

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

extern "C" {
  void zgbsv_ ( const int& N, const int& KL, const int& KU, const int& NRHS,
		Complex** AB, const int& LDAB, int* IPIV,
		Complex**  B, const int& LDB, int& INFO );
};

// 複素 3重 帯称行列 A の線形方程式 A x = b を解く簡易関数
// Input: A[3][N], b[N]  Output: x[N]
// A[1][N] に対角成分、A[0][N] に下副対角成分、A[2][N] に上副対角成分
//
template <int N> int zgbsv3( Complex A[3][N], Complex x[N], Complex b[N] )
{
  int i, j, info;
  static int ipiv[N];
  static Complex AB[N][N];

  for( j=0; j<N; j++ ){
    for( i=0; i<3; i++ ){
      AB[j][1+i] = A[i][j];
    }
    x[j] = b[j];
  }

  zgbsv_( N, 1, 1, 1, (Complex**)AB, N, ipiv, (Complex**)x, N, info );

  return info;
}


// 実 5重 帯称行列 A の線形方程式 A x = b を解く簡易関数
// Input: A[3][N], b[N]  Output: x[N]
// A[2][N] に対角成分、A[0,1][N] に上副対角成分、A[3,4][N] に下副対角成分
//
template <int N> int zgbsv5( Complex A[5][N], Complex x[N], Complex b[N] )
{
  int i, j, info;
  static int ipiv[N];
  static Complex AB[N][N];

  for( j=0; j<N; j++ ){
    for( i=0; i<5; i++ ){
      AB[j][2+i] = A[i][j];
    }
    x[j] = b[j];
  }

  zgbsv_( N, 2, 2, 1, (Complex**)AB, N, ipiv, (Complex**)x, N, info );

  return info;
}


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

  Complex A[5][N];
  Complex x[N], b[N];

  for( j=0; j<N; j++ ){
    A[0][j] = Complex(1.0, 2.0);
    A[1][j] = Complex(1.0, 2.0);
    A[2][j] = Complex(-2.0, 0.0);
    A[3][j] = Complex(1.0, 2.0);
    A[4][j] = Complex(1.0, 2.0);
       b[j] = Complex(4.0, 2.0);
  }

  info = zgbsv5( A, x, b );

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

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

  printf("# Error.\n");
  for( i=0; i<N; i++ ){
    Complex sum=0.0;
    for( j=i-2; j<=i+2; j++ ){
      if( j<0 || N<=j ) continue;
      sum += A[2+i-j][j]*x[j];
    }
    sum -= b[i];

    printf("%+f%+f\n", real(sum), imag(sum) );
  }

  return 0;
}