// dgtsv.cc
// 実3重対角行列の線形方程式を解く

#include <stdio.h>

extern "C" {
  void dgtsv_ ( const int& N, const int& NRHS,
		double* DL, double* D, double* DU,
		double** 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 dgtsv( double A[3][N], double x[N], double b[N] )
{
  int j, info;
  static double DL[N], D[N], DU[N];

  for( j=0; j<N; j++ ){
    DU[j] = A[0][j+1];
     D[j] = A[1][j];
    DL[j] = A[2][j];

    x[j] = b[j];
  }

  dgtsv_( N, 1, DL, D, DU, (double**)x, N, info );

  return info;
}


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

  double A[3][N];
  double x[N], b[N];

  for( j=0; j<N; j++ ){
    A[0][j] = 1.0;
    A[1][j] = -2.0;
    A[2][j] = 1.0;
       b[j] = 5.0;
  }

  info = dgtsv( A, x, b );

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

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

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

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

  return 0;
}