/*
 * cg.C --- 入出力に iostream を使う
 *    g++ -O -o cg cg.C  -L/usr/local/lib -lmatrix
 */

#include <iostream.h>
#include <cmath>
extern "C" {
#include <matrix.h>
}

// m行n列の零行列を作って返す
matrix zeros(int m, int n)
{
  int i, j;
  matrix tmp;
  tmp = new_matrix(m, n);
  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++)
      tmp[i][j] = 0.0;
  return tmp;
}

// m行n列の行列 z を零行列にする
void zeros(int m, int n, matrix z)
{
  int i, j;
  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++)
      z[i][j] = 0.0;
}

// n次元の零ベクトルを作って返す
vector zeros(int n)
{
  int i;
  vector tmp;
  tmp = new_vector(n);
  for (i = 0; i < n; i++)
    tmp[i] = 0.0;
  return tmp;
}

// n次元のベクトル z を零ベクトルにする
void zeros(int n, vector z)
{
  int i;
  for (i = 0; i < n; i++)
    z[i] = 0.0;
}

// n 次元ベクトルの内容を表示する
void print_vector(int n, vector x)
{
  int i;
  for (i = 0; i < n; i++) {
    cout << " " << x[i];
    if (i % 5 == 4)
      cout << endl;
  }
  if (i % 5 != 0)
    cout << endl;
}

// m行n列の行列の内容を表示する
void print_matrix(int m, int n, matrix a)
{
  int i, j;
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      cout << a[i][j];
      if (j % 5 == 4)
        cout << endl;
    }
    if (j % 5 != 0)
      cout << endl;
  }
}

// ベクトル x, y の内積 (x,y) を計算する
double dotproduct(int n, vector x, vector y)
{
  int i;
  double sum;
  sum = 0;
  for (i = 0; i < n; i++)
    sum += x[i] * y[i];
  return sum;
}

// ベクトル x のノルム ||x|| を計算する
double norm(int n, vector x)
{
  return sqrt(dotproduct(n, x, x));
}

// 行列 A, ベクトル b の積 A b を計算して、ベクトル Ab に書き込む
void multiply_mv(int m, int n, vector Ab, matrix A, vector b)
{
  int i;
  for (i = 0; i < m; i++)
      Ab[i] = dotproduct(n, A[i], b);
}

// ベクトル x をベクトル y にコピーする
void copy_vector(int n, vector y, vector x)
{
  int i;
  for (i = 0; i < n; i++)
    y[i] = x[i];
}

// CG 法で n 元連立1次方程式 A x=b を解く。
// 呼び出すときの x は初期ベクトル, 戻ってきたときの x は解となる。
// ||残差|| < ε ||b|| となったら反復を停止する。
// 反復回数を maxiter に返す。
void cg(int n, matrix A, vector b, vector x, double eps, int *maxiter)
{
  int i, k;
  double eps_b, pAp, alpha, beta;
  vector r, p, Ap;

  r = new_vector(n);
  p = new_vector(n);
  Ap = new_vector(n);
  /* eps_b := ε ||b|| */
  eps_b = eps * norm(n, b);
  /* r := A x */
  multiply_mv(n, n, r, A, x);
  /* r := b - A x */
  for (i = 0; i < n; i++)
    r[i] = b[i] - r[i];
  /* p := r */
  copy_vector(n, p, r);
  /* 反復 */
  for (k = 0; k < n; k++) {
    /* Ap := A * p */
    multiply_mv(n, n, Ap, A, p);
    /* pAp := (p, Ap) */
    pAp = dotproduct(n, p, Ap);
    /* α = (r,p)/(p,Ap) */
    alpha = dotproduct(n, r, p) / pAp;
    /* x, r の更新 */
    for (i = 0; i < n; i++) {
      /* x = x + α*p */
      x[i] += alpha * p[i];
      /* r = r - α*A*p */
      r[i] -= alpha * Ap[i];
    }
    /* ||r||<ε||b|| ならば反復を終了 */
    if (norm(n, r) < eps_b)
      break;
    /* β := - (r, Ap) / (p, Ap) */
    beta = - dotproduct(n, r, Ap) / pAp;
    /* p := r + β*p */
    for (i = 0; i < n; i++)
      p[i] = r[i] + beta * p[i];
  }
  if (maxiter != NULL)
    *maxiter = k;
}

int main()
{
  int n, N, i, niter;
  matrix a;
  vector x, b;
  double c1, c2;
  /* 問題のサイズを決める */
  cout << "c1, c2, N=";
  cin >> c1 >> c2 >> N;
  n = N / 2;
  N = 2 * n;
  /* 行列 A, ベクトル x, b を記憶する変数の準備 */
  a = zeros(N, N);
  x = new_vector(N);
  b = new_vector(N);
  /* 係数行列を決める */
  for (i = 0; i < N; i++) {
    a[i][i] = 4.0;
    if (i != 0)
      a[i][i-1] = -1.0;
    if (i != N-1)
      a[i][i+1] = -1.0;
  }
  a[n-1][n] = 0.0;
  a[n][n-1] = 0.0;
  /* */
  for (i = 0; i < n; i++) {
    if (i != 0)
      a[i][i-1] *= c1;
    a[i][i] *= c1;
    a[i][i+1] *= c1;
  }
  for (i = n; i < N; i++) {
    a[i][i-1] *= c2;
    a[i][i] *= c2;
    if (i != N-1)
      a[i][i+1] *= c2;
  }
  /* 係数行列を表示する */
  if (N < 10) {
    cout << "a=" << endl;
    print_matrix(N, N, a);
  }
  /* 解 x を決める */
  for (i = 0; i < N; i++)
    x[i] = i;
  /* 解に対応する右辺 b を計算する */
  multiply_mv(N, N, b, a, x);
  /* b を表示する */
  if (N < 10) {
    cout << "b=" << endl;
    print_vector(N, b);
  }
  /* 初期ベクトル x を零ベクトルにする */
  zeros(N, x);
  /* CG 法で解く */
  cg(N, a, b, x, 1e-12, &niter);
  /* 解を表示する */
  cout << "x=" << endl;
  print_vector(N, x);
  cout << "反復回数=" << niter << endl;

  return 0;
}